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ABSTRACT 

We test the statistical isotropy and Gaussianity of the cosmic microwave background (CMB) anisotropies using ob¬ 
servations made by the Planck satellite. Our results are based mainly on the full Planck mission for temperature, 
but also include some polarization measurements. In particular, we consider the CMB anisotropy maps derived from 
the multi-frequency Planck data by several component-separation methods. For the temperature anisotropies, we find 
excellent agreement between results based on these sky maps over both a very large fraction of the sky and a broad 
range of angular scales, establishing that potential foreground residuals do not affect our studies. Tests of skewness, 
kurtosis, multi-normality, iV-point functions, and Minkowski functionals indicate consistency with Gaussianity, while 
a power deficit at large angular scales is manifested in several ways, for example low map variance. The results of a 
peak statistics analysis are consistent with the expectations of a Gaussian random field. The “Cold Spot” is detected 
with several methods, including map kurtosis, peak statistics, and mean temperature profile. We thoroughly probe the 
large-scale dipolar power asymmetry, detecting it with several independent tests, and address the subject of a poste¬ 
riori correction. Tests of directionality suggest the presence of angular clustering from large to small scales, but at a 
significance that is dependent on the details of the approach. We perform the first examination of polarization data, 
finding the morphology of stacked peaks to be consistent with the expectations of statistically isotropic simulations. 
Where they overlap, these results are consistent with the Planck 2013 analysis based on the nominal mission data and 
provide our most thorough view of the statistics of the CMB fluctuations to date. 

Key words, cosmology: observations - cosmic background radiation - polarization - methods: data analysis - methods: 
statistical 
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1. Introduction 

This paper, one of a set associated with the 2015 release 
of data from the Planck * 1 mission (Planck Collaboration I 
2015), describes a set of studies undertaken to determine 
the statistical properties of both the temperature and po¬ 
larization anisotropies of the cosmic microwave background 
(CMB). 

The standard cosmological model is described well by 
the Friedmann-Lemaitre-Robertson-Walker solution of the 
Einstein field equations. This model is characterized by a 
homogeneous and isotropic background metric and a scale 
factor of the expanding Universe. It is hypothesized that 
at very early times the Universe went through a period 
of accelerated expansion, the so-called “cosmological infla¬ 
tion,” driven by a hypothetical scalar field, the “inflaton.” 
During inflation the Universe behaves approximately as a 
de Sitter space, providing the conditions by which some of 
its present properties can be realized and specifically re¬ 
laxing the problem of initial conditions. In particular, the 
seeds that gave rise to the present large-scale matter distri¬ 
bution via gravitational instability originated as quantum 
fluctuations of the inflaton about its vacuum state. These 
fluctuations in the inflaton produce energy density pertur¬ 
bations that are distributed as a statistically homogeneous 
and isotropic Gaussian random field. Linear theory relates 
those perturbations to the temperature and polarization 
anisotropies of the CMB, implying a distribution for the 
anisotropies very close to that of a statistically isotropic 
Gaussian random field. 

The aim of this paper is to use the full mission Planck 
data to test the Gaussianity and isotropy of the CMB as 
measured in both intensity and, in a more limited capacity, 
polarization. Testing these fundamental properties is cru¬ 
cial for the validation of the standard cosmological scenario, 
and has profound implications for our understanding of the 
physical nature of the Universe and the initial conditions 
of structure formation. Moreover, the confirmation of the 
statistically isotropic and Gaussian nature of the CMB is 
essential for justifying the corresponding assumptions usu¬ 
ally made when estimating the CMB power spectra and 
other quantities to be obtained from the Planck data. In¬ 
deed, the isotropy and Gaussianity of the CMB anisotropies 
are implicitly assumed in critical science papers from the 
2015 release, in particular those describing the likelihood 
and the derivation of cosmological parameter constraints 
(Planck Collaboration XI 2015; Planck Collaboration XIII 
2015). Conversely, if the detection of significant deviations 
from these assumptions cannot be traced to known system¬ 
atic effects or foreground residuals, the presence of which 
should be diagnosed by the statistical tests set forth in 
this paper, this would necessitate a major revision of the 
current methodological approaches adopted in deriving the 
mission’s many science results. 


* Corresponding author: A. J. Banday anthony. banday@irap. 
omp.eu 

1 Planck (http://www.esa.int/Planck) is a project of the Eu¬ 
ropean Space Agency (ESA) with instruments provided by two 
scientific consortia funded by ESA member states and led by 
Principal Investigators from France and Italy, telescope reflec¬ 
tors provided through a collaboration between ESA and a sci¬ 
entific consortium led and funded by Denmark, and additional 
contributions from NASA (USA). 


Well-understood physical processes due to the inte¬ 
grated Sachs-Wolfe (ISW) effect (Planck Collaboration 
XVII 2014; Planck Collaboration XXI 2015) and gravita¬ 
tional lensing (Planck Collaboration XIX 2014; Planck Col¬ 
laboration XV 2015) lead to secondary anisotropies that 
exhibit marked deviation from Gaussianity. In addition, 
Doppler boosting, due to our motion with respect to the 
CMB rest frame, induces both a dipolar modulation of 
the temperature anisotropies and an aberration that cor¬ 
responds to a change in the apparent arrival directions of 
the CMB photons (Challinor & van Leeuwen 2002). Both 
of these effects are aligned with the CMB dipole, and were 
detected at a statistically significant level on small angular 
scales in Planck Collaboration XXVII (2014). Beyond these, 
Planck Collaboration XXIII (2014, hereafter PCIS13) es¬ 
tablished that the Planck 2013 data set showed little evi¬ 
dence for non-Gaussianity, with the exception of a number 
of CMB temperature anisotropy anomalies on large angu¬ 
lar scales that confirmed earlier claims based on WMAP 
data. Moreover, given that the broader frequency cover¬ 
age of the Planck instruments allowed improved compo¬ 
nent separation methods to be applied in the derivation of 
foreground-cleaned CMB maps, it was generally considered 
that the case for anomalous features in the CMB had been 
strengthened. Hence, such anomalies have attracted consid¬ 
erable attention in the community, since they could be the 
visible traces of fundamental physical processes occurring 
in the early Universe. 

However, the literature also supports an ongoing debate 
about the significance of these anomalies. The central issue 
in this discussion is connected with the role of a posteri¬ 
ori choices — whether interesting features in the data bias 
the choice of statistical tests, or if arbitrary choices in the 
subsequent data analysis enhance the significance of the fea¬ 
tures. Indeed, the WMAP team (Bennett et al. 2011) base 
their rejection of the presence of anomalies in the CMB on 
such arguments. Of course, one should attempt to correct 
for any choices that were made in the process of detect¬ 
ing an anomaly. However, in the absence of an alternative 
model for comparison to the standard Gaussian, statisti¬ 
cally isotropic one adopted to quantify significance, this is 
often simply not possible. In this work, whilst it is recog¬ 
nized that care must be taken in the assessment of signif¬ 
icance, we proceed on the basis that allowing a posteriori 
reasoning permits us to challenge the limits of our existing 
knowledge (Pontzen & Peiris 2010). That is, by focusing 
on specific properties of the observed data that are shown 
to be empirically interesting, we may open up new paths 
to a better theoretical understanding of the Universe. We 
will clearly describe the methodology applied to the data, 
and attempt to study possible links among the anomalies 
in order to search for a physical interpretation. 

The analysis of polarization data introduces a new op¬ 
portunity to explore the statistical properties of the CMB 
sky, including the possibility of improvement of the sig¬ 
nificance of detection of large-scale anomalies. However, 
this cannot be fully included in the current data assess¬ 
ment, since the component-separation products in polar¬ 
ization are high-pass filtered to remove large angular scales 
(Planck Collaboration IX 2015), owing to the persistence of 
significant systematic artefacts originating in the High Fre¬ 
quency Instrument (HFI) data (Planck Collaboration VII 
2015; Planck Collaboration VIII 2015). In addition, limi¬ 
tations of the simulations with which the data are to be 


Article number, page 2 of 66 



Planck Collaboration: Isotropy and statistics of the CMB 


compared (Planck Collaboration XII 2015), in particular a 
significant mismatch in noise properties, limit the extent to 
which any polarization results can be included. Therefore, 
we only present a stacking analysis of the polarized data, al¬ 
though this is a significant extension of previous approaches 
found in the literature. 

With future Planck data releases, it will be important 
to determine in more detail whether there are any pecular- 
ities in the CMB polarization, and if so, whether they are 
related to existing features in the CMB temperature field. 
Conversely, the absence of any corresponding features in 
polarization might imply that the the temperature anoma¬ 
lies (if they are not simply flukes) could be due to a sec¬ 
ondary effect such as the ISW effect, or alternative scenarios 
in which the anomalies arise from physical processes that 
do not correlate with the temperature, e.g., texture or de¬ 
fect models. Either one of these possible outcomes could 
yield a breakthrough in understanding the nature of the 
CMB anomalies. Of course, there also remains the possibil¬ 
ity that anomalies may be found in the polarization data 
that are unrelated to existing features in the temperature 
measurements. 

Following the approach established in Planck Collabo¬ 
ration XXIII (2014), throughout this paper we quantify the 
significance of a test statistic in terms of the p- value. This 
is the probability of obtaining a test statistic at least as ex¬ 
treme as the observed one, under the assumption that the 
null hypothesis (i.e., primordial Gaussianity and isotropy of 
the CMB) is true. In some tests, where it is clearly justified 
to only use a one-tailed probability, the p -value is replaced 
by the corresponding upper- or lower-tail probability. 

This paper covers all relevant aspects related to the phe¬ 
nomenological study of the statistical isotropy and Gaus¬ 
sian nature of the CMB measured by the Planck satel¬ 
lite. Specific theoretically-motivated model constraints on 
isotropy or non-Gaussianity, as might arise from non¬ 
standard inflationary models, the geometry and topology 
of the Universe, and primordial magnetic fields are pro¬ 
vided in the companion papers (Planck Collaboration XVII 
2015; Planck Collaboration XX 2015; Planck Collaboration 
XVIII 2015; Planck Collaboration XIX 2015). The paper is 
organized as follows. Section 2 summarizes the Planck full 
mission data used for the analyses, and important limita¬ 
tions of the polarization maps that are studied. Section 3 
describes the characteristics of the simulations that consti¬ 
tute our reference set of Gaussian sky maps representative 
of the null hypothesis. In Sect. 4 the null hypothesis is tested 
with a number of standard tests that probe different aspects 
of non-Gaussianity. Several important anomalous features 
of the CMB sky, originally detected with the WMAP data 
and subsequently confirmed in PCIS13, are reassessed in 
Sect. 5. Aspects of the CMB fluctuations specifically re¬ 
lated to dipolar asymmetry are examined in Sect. 6. The 
sensitivity of the results for a number of statistical tests to 
the sky fraction is examined in Sect. 7. Section 8 presents 
tests of the statistical nature of the polarization signal ob¬ 
served by Planck using a local analysis of stacked patches 
of the sky. Finally, Sect. 9 provides the main conclusions of 
the paper. 

2. Data description 

In this paper, we use data from the Planck-2015 full mis¬ 
sion data release. This contains approximately 29 months 


of data for the HFI and 50 months for the Low Frequency 
Instrument (LFI). The release includes sky maps at nine 
frequencies in intensity (seven in polarization), with corre¬ 
sponding “half-mission” maps that are generated by split¬ 
ting the full-mission data sets in various ways. The maps 
are provided in HEALPix format (Gorski et al. 2005), 2 with 
a pixel size defined by the V si d e parameter. This set of maps 
allows a variety of consistency checks to be made, together 
with estimates of the instrumental noise contributions to 
our analyses and limits on time-varying systematic arte¬ 
facts. Full details are provided in a series of companion 
papers (Planck Collaboration II 2015; Planck Collabora¬ 
tion III 2015; Planck Collaboration IV 2015; Planck Col¬ 
laboration V 2015; Planck Collaboration VI 2015; Planck 
Collaboration VII 2015; Planck Collaboration VIII 2015). 

Our main results are based on estimates of the CMB 
generated by four distinct component-separation algo¬ 
rithms — Commander, NILC, SEVEM, and SMICA — as de¬ 
scribed in Planck Collaboration IX (2015). These effectively 
combine the raw Planck frequency maps in such a way as to 
minimize foreground residuals from diffuse Galactic emis¬ 
sion. Note that the additional information in the full mis¬ 
sion data set allows us to improve the reconstruction noise 
levels by roughly a factor of 2 (in temperature) as com¬ 
pared to the Planck -2013 nominal mission data release. The 
CMB intensity maps were derived using all channels, from 
30 to 857 GHz, and provided at a common angular resolu¬ 
tion of 5' FWHM and N s ^ e = 2048. The intensity maps 
are only partially corrected for the second order tempera¬ 
ture quadrupole (Kamionkowski & Knox 2003). Therefore, 
where appropriate, the component-separated maps should 
be corrected for the residual contribution (Notari & Quar- 
tin 2015), specifically as described in Planck Collaboration 
IX (2015). The polarization solutions include all channels 
sensitive to polarization, from 30 to 353 GHz, at a resolution 
of 10' FWHM and V si d e = 1024. Possible residual emission 
is then mitigated in our analyses by the use of sky-coverage 
masks, provided for both intensity and polarization. 

Since in some cases it is important to study the fre¬ 
quency dependence of the cosmological signal, either to es¬ 
tablish its primordial origin or to test for the frequency 
dependence associated with specific effects such as Doppler 
boosting (see Sect. 6.4), we also consider the foreground- 
cleaned versions of the 100, 143, and 217 GHz sky maps 
generated by the SEVEM algorithm (Planck Collaboration 
IX 2015), hereafter referred to as SEVEM-100, SEVEM-143, 
and SEVEM-217, respectively. 

For the present release, a post-processing high-pass¬ 
filtering has been applied to the CMB polarization maps 
in order to mitigate residual large-scale systematic errors 
in the HFI channels (Planck Collaboration VII 2015). The 
filter results in the elimination of structure in the maps on 
angular scales larger than about 10°, and a weighted sup¬ 
pression of power down to scales of 5°, below which the 
maps remain unprocessed. 

Lower-resolution versions of these data sets are also used 
in the analyses presented in this paper. The downgrading 
procedure for maps is to decompose them into spherical har¬ 
monics on the full sky at the input HEALPix resolution. The 
spherical harmonic coefficients, a£ m , are then convolved to 


2 http://healpix.sourceforge.net 
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Table 1 . Standardized data sets used in this paper. The resolu¬ 
tions of the sky maps used are defined in terms of the N S id e pa¬ 
rameter and corresponding FWHM of the Gaussian beam with 
which they are convolved. The corresponding common masks 
and the fraction of unmasked pixels used for analysis are also 
specified. 


FWHM Mask Unmasked 


Wide [arcmin] pixels [%] 


2048 . 5 UT78 77.6 

2048 . 5 UTA76 76.1 

1024 . 10 UTio24 7 6 75.6 

512. 20 UT 5 i 2 74 73.7 

256 . 40 UT 256 73 72.5 

128 . 80 UTi 28 70 69.7 

64 . 160 UT 64 6 7 67.0 

32 . 320 UT 32 6 4 63.8 

16 . 640 UTi 6 58 58.4 

1024 . 10 UPB77 77.4 


the new resolution using 


a 


out 

im 


b out p out 

w' 


Um ) 


(i) 


where bn is the beam transfer function, p£ is the HEALPix 
pixel window function, and the “in” and “out” superscripts 
denote the input and output resolutions. They are then 
synthesized into a map directly at the output HEALPix res¬ 
olution. Masks are downgraded in a similar way. The bi¬ 
nary mask at the starting resolution is first downgraded 
like a temperature map. The smooth downgraded mask is 
then thresholded by setting pixels where the value is less 
than 0.9 to zero and all others to unity in order to make 
a binary mask. Table 1 lists the Wide and FWHM values 
defining the resolution of these maps, together with the dif¬ 
ferent masks and their sky coverages that accompany the 
signal maps. In general, we make use of standardized masks 
that are the union of those associated with the individual 
component-separation methods. 

As recommended in Planck Collaboration IX (2015), 
the mask UT78 is adopted for all high-resolution analyses 
of temperature data. UTA76 is an extended version of this 
mask more suitable for some non-Gaussianity studies. The 
mask preferred for polarization studies, UPB77, is again the 
union of those determined for each component separation 
method, but in addition the polarized point sources de¬ 
tected at each frequency channel are excluded. These masks 
are then downgraded for lower-resolution studies. As a con¬ 
sequence of the common scheme applied in order to gen¬ 
erate such low-resolution masks, they are generally more 
conservative than the corresponding ones used in the 2013 
analyses. 

In what follows, we will undertake analyses of the data 
at a given resolution denoted by a specific WM e va l ue - Un¬ 
less otherwise stated, this implies that the data have been 
smoothed to a corresponding FWHM as described above, 
and a standardized mask employed. Irrespective of the res¬ 
olution in question, we will then often simply refer to the 
latter as the “common mask.” 


3. Simulations 

The results presented in this paper are derived using the 
extensive full focal plane (FFP8) simulations described in 


Planck Collaboration XII (2015). Of most importance are 
the Monte Carlo (MC) simulations that provide the refer¬ 
ence set of Gaussian sky maps used for the null tests em¬ 
ployed here. They also form the basis of any debiasing in the 
analysis of the real data as required by certain statistical 
methods. 

The simulations include both CMB signal and instru¬ 
mental noise realizations that capture important character¬ 
istics of the Planck scanning strategy, telescope, detector 
responses, and data reduction pipeline over the full mis¬ 
sion period. In particular, the signal realizations include 
FEBeCoP (Mitra et al. 2011) beam convolution at each of the 
Planck frequencies, and are propagated through the various 
component-separation pipelines using the same weights as 
derived from the Planck full mission data analysis. 

The FFP8 fiducial CMB power spectrum has been 
adopted from our best estimate of the cosmological pa¬ 
rameters from the first Planck data release (Planck Col¬ 
laboration I 2014). This corresponds to a cosmology with 
baryon density given by c^b = = 0.0222, cold dark 

matter density uj c = f} c h 2 = 0.1203, neutrino energy den¬ 
sity = fl u h 2 = 0.00064, density parameter for the 
cosmological constant = 0.6823, Hubble parameter 
Ho = 100/ikms -1 Mpc -1 with h = 0.6712, spectral index 
of the power spectrum of the primordial curvature pertur¬ 
bation n s = 0.96, and amplitude of the primordial power 
spectrum (at k = 0.05 Mpc -1 ) A s = 2.09 x 10 -9 , and with 
the Thomson optical depth through reionization defined to 
be r = 0.065. Each realization of the CMB sky is generated 
including lensing, Rayleigh scattering, and Doppler boost¬ 
ing effects, the latter two of which are frequency-dependent. 
A second order temperature quadrupole (Kamionkowski & 
Knox 2003) is added to each simulation with an ampli¬ 
tude corresponding to the residual uncorrected contribu¬ 
tion present in the observed data, as described in Planck 
Collaboration XII (2015). 

However, the Planck maps were effectively renormalized 
by approximately 2 % to 3 % in power in the time between 
the generation of the FFP8 simulations and the final maps. 
As discussed in Planck Collaboration XII (2015), correction 
for this calibration effect should have no significant impact 
on cosmological parameters. As recommended, in this paper 
the CMB component of the simulations is simply rescaled 
by a factor of 1.0134 before analysis. 

Of somewhat more importance is an observed noise 
mismatch between the simulations and the data. Whilst 
this has essentially no impact on studies of temperature 
anisotropy, it imposes important limitations on the statis¬ 
tical studies of polarization sky maps that can be included 
here. Conversely, analyses based on 1-point statistics, such 
as the variance, and the TV-point correlation functions have 
played important roles in establishing the nature of this 
mismatch, which seems to be scale-dependent with an am¬ 
plitude around 20 % at lower resolutions but falling to a 
few per cent at higher resolution. As a consequence, this 
paper only includes results from a stacking analysis of the 
polarized data, in which the stacking of the data themselves 
necessarily acts to lower the effect of the noise mismatch. 
Polarization studies that do not rely on auto-statistics can 
still yield interesting new results, as found in Planck Col¬ 
laboration XIII (2015); Planck Collaboration XVII (2015); 
Planck Collaboration XVIII (2015). 
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4. Tests of non-Gaussianity 

There is no unique signature of non-Gaussianity, but the 
application of a variety of tests over a range of angular 
scales allows us to probe the data for departures from the¬ 
oretically motivated Gaussian statistics. One of the more 
important tests in the context of inflationary cosmology 
is related to the analysis of the bispectrum. This is dis¬ 
cussed thoroughly in Planck Collaboration XVII (2015), 
and is therefore not discussed further in this paper. In this 
section, we present the results from a variety of statistical 
tools. Unless otherwise specified, the analyses are applied to 
all four component separation products (Commander, NILC, 
SEVEM, and SMICA) at a given resolution with the accompa¬ 
nying common mask, and significance levels are determined 
by comparison with the corresponding results derived from 
the FFP8 simulations, with typically 1000 being used for 
this purpose. Establishing the consistency of the results de¬ 
rived from the different component-separation techniques is 
essential in order to be able to make robust claims about the 
statistical nature of the observed temperature fluctuations, 
and potential deviations from Gaussianity. 
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Fig. 1. Variance, skewness, and kurtosis for the four different 
component-separation methods — Commander (red), NILC (or¬ 
ange), SEVEM (green), and SMICA (blue) — compared to the dis¬ 
tributions derived from 1000 Monte Carlo simulations. 


4.1. One-dimensional moments 


In this section we consider simple tests of Gaussianity based 
on the variance, skewness, and kurtosis of the CMB tem¬ 
perature maps. Previous analyses found an anomalously low 
variance in the WMAP sky maps (Monteserfn et al. 2008; 
Cruz et al. 2011), which was subsequently confirmed in an 
analysis of the Planck 2013 data (PCIS13). 

Cruz et al. (2011) developed the unit variance estimator 
to determine the variance, <7 q, of the CMB signal on the sky 
in the presence of noise. The normalized CMB map, u x , is 
given by 


u f ( a x,o) 





( 2 ) 


where Xi is the observed temperature at pixel i and of noise 
is the noise variance for that pixel. Although this esti¬ 
mator is not optimal, Cruz et al. (2011) and Monteserfn 
et al. (2008) have demonstrated that it is unbiased and suf¬ 
ficiently accurate for our purposes. The noise variance is 
estimated from the noise simulations for each component- 
separation algorithm. The CMB variance is then estimated 
by requiring that the variance of the normalized map u x 
is unity. The skewness and kurtosis can then be obtained 
from the appropriately normalized map. 

Figure 1 presents results for the variance, skewness, and 
kurtosis determined from the data at a resolution of 5', 
Aside = 2048. Good agreement between the component sep¬ 
aration techniques is found, with small discrepancies likely 
due to sensitivity to the noise properties and their variation 
between methods. 

Table 2 summarizes the lower-tail probabilities, defined 
as the percentage of MC simulations that show a lower vari¬ 
ance, skewness, or kurtosis than the observed map, for these 
analyses. The results are in good agreement with PCIS13; 
the skewness and kurtosis are compatible with simulations, 
but the variance is marginally lower than in the simulations. 

Although the variance is observed to be low, the re¬ 
sults could still be affected by the presence of residual fore¬ 
grounds at small scales in these maps, so that the true vari- 


Table 2. Lower-tail probabilities for the variance, skewness, and 
kurtosis of the component-separated maps. 


Probability [%] 


Method Variance Skewness Kurtosis 


Commander .... 

3.2 

17.2 

35.3 

NILC. 

3.3 

20.9 

30.9 

SEVEM . 

1.9 

20.5 

56.8 

SMICA . 

1.4 

21.1 

48.2 

SEVEM-100 .... 

3.4 

13.4 

67.5 

SEVEM-143 _ 

2.4 

16.9 

61.2 

SEVEM-217 .... 

3.4 

11.4 

58.3 


Table 3. Lower-tail probabilities for the iV-pdf y 2 statistics 
derived from the Planck 2015 component-separated maps at 
A S ide = 16 and 32. 


Probability [%] 


Aside Comm. NILC SEVEM SMICA 

16 . 24.7 26.2 25.4 24.5 

32 . 11.9 20.8 10.6 10.8 


ance would be lower still. We assess this by application of 
the estimator to the cleaned frequency maps SEVEM-100, 
SEVEM-143, and SEVEM-217. The results, also presented in 
Table 2, are similar to those found for the combined map, 
although slightly less significant, which is most likely at¬ 
tributable to higher noise in the cleaned frequency maps. 

In conclusion, a simple statistical assessment of the 
Planck 2015 data using skewness and kurtosis shows no 
evidence for non-Gaussianity, although a low variance is 
found, which we will readdress in Sect. 5.1. 

4.2. Testing the multi-normality of the CMB 

Under the assumption of Gaussianity, the probability den¬ 
sity function (PDF) of the TV-dimensional pixelized tem¬ 
perature map is given by a multivariate Gaussian function: 
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f(T) = 


(27r) JVpix/2 det CV2 


exp 


4( rc ' 


i T T j 


(3) 


where T is a vector formed from the measured temperatures 
T{x) over all positions allowed by the applied mask, A pix is 
the number of pixels in the vector, and C is the covariance 
of the Gaussian field (of size A pix x A pix ). 

Although the calculation of TC -1 T t can be achieved by 
conjugate gradient methods, the evaluation of det C remains 
computationally difficult for the full Planck resolution at 
HEALPix A^ide = 2048. At a lower resolution, the problem 
is tractable, and the noise level can also be considered negli¬ 
gible compared to the CMB signal. That implies that under 
the assumption of isotropy the covariance matrix C is fully 
defined by the Planck angular power spectrum (Ct): 


C 


1=2 


2^ + 1 
47r 


CebjPe (cosOij ), 


(4) 


where Qj is the covariance between pixels i and j, Oij is 
the angle between them, Pg are the Legendre polynomials, 
hi is an effective window function describing the combined 
effects of the instrumental beam and pixel window at reso¬ 
lution A S ide, and ^ max is the maximum multipole probed. 

Under the multivariate Gaussian hypothesis, the argu¬ 
ment of the exponential in Eq. (3) should follow a y 2 distri¬ 
bution with A pix degrees of freedom, or, equivalently (for 
A P i x 1) a normal distribution Af (A pix , Y / 2A pix ). 

These y 2 statistics are computed for the Planck 2015 
component-separated CMB maps at A si d e = 16 and 32, 
then compared with the equivalent quantities derived from 
the corresponding FFP8 simulations. For those cases in 
which the covariance matrix is ill-conditioned, we use a 
principal component analysis approach to remove the low¬ 
est degenerate eigenvalues of the covariance matrix (see, 
e.g., Curto et al. 2011). This process is equivalent to adding 
uncorrelated regularization noise of amplitude « 1 fiK to 
the data before inversion. The results of the analysis are 
presented in Table 3 and indicate that the data are consis¬ 
tent with Gaussianity. We note that the lower-tail probabil¬ 
ities for the A-pdf decrease when the resolution of the data 
is increased from A si d e = 16 to 32. However, this behaviour 
is consistent with that seen for simulations, and should not 
be considered to be significant. 


4.3. N-point correlation functions 

In this section, we present tests of the non-Gaussianity of 
the Planck 2015 temperature CMB maps using real-space 
A-point correlation functions. While harmonic-space meth¬ 
ods are often preferred over real-space methods for study¬ 
ing primordial fluctuations, real-space methods have an ad¬ 
vantage with respect to systematic errors and foregrounds, 
since such effects are usually localized in real space. It is 
therefore important to analyse the data in both spaces in 
order to highlight different features. 

An A-point correlation function is defined as the aver¬ 
age product of A temperatures, measured in a fixed relative 
orientation on the sky, 


(5) 


where the unit vectors ni,..., tin span an A-point poly¬ 
gon. Under the assumption of statistical isotropy, these 
functions depend only on the shape and size of the A-point 
polygon, and not on its particular position or orientation 
on the sky. Hence, the smallest number of parameters that 
uniquely determines the shape and size of the A-point poly¬ 
gon is 2A — 3. 

The correlation functions are estimated by simple prod¬ 
uct averages over all sets of A pixels fulfilling the geometric 
requirements set by 0i,..., 0 2 n-3 characterizing the shape 
and size of the polygon, 


C N (0 1, • • • ,#2iV— 3 ) = 


Ei «■■■<) pi• • • 


Ei 


W\ 


■ W 


(6) 


N) 


Pixel weights w \,.. *, w l N can be introduced in order to re¬ 
duce noise or mask boundary effects. Here they represent 
masking by being set to 1 for included pixels and to 0 for 
excluded pixels. 

The shapes of the polygons selected for the analysis are 
the pseudo-collapsed and equilateral configurations for the 

3- point function, and the rhombic configuration for the 4- 
point function, composed of two equilateral triangles that 
share a common side. We use the same definition of pseudo- 
collapsed as in Eriksen et al. (2005), i.e., an isosceles trian¬ 
gle where the length of the baseline falls within the second 
bin of the separation angles. The length of the longer edge 
of the triangle, #, parameterizes its size. Analogously, in the 
case of the equilateral triangle and rhombus, the size of the 
polygon is parameterized by the length of the edge, 0. Note 
that these functions are chosen for ease of implementation, 
not because they are better suited for testing Gaussianity 
than other configurations. For a Gaussian field, Wick’s the¬ 
orem (Wick 1950) means that the ensemble average of the 

4- point function may be written in terms of the 2-point 
function. In the following, all results refer to the connected 
4-point function, i.e., are corrected for this Gaussian con¬ 
tribution. 

We use a simple y 2 statistic to quantify the agreement 
between the observed data and simulations, defined by 


iVbin 

X 2 = F (pN(0i) - (Cjv(0»))) (pN(0j) - (CnV,))) ■ 

fj =1 

(7) 

Here, Cn{0i) is the A-point correlation function for the 
bin with separation angle 6^, (CV(#i)) is the correspond¬ 
ing average from the MC simulation ensemble, and Abi n is 
the number of bins used for the analysis. If C%{0i) is the 
fcth simulated A-point correlation function and A sim is the 
number of simulations, then the covariance matrix M^- is 
given by 


1 % = w~ E (A fe) (P - <cw(0o>) x 

sim k=1 

(c^\ej)-{c jv (%)>) , ( 8 ) 

where A s ' im = A sim — 1. Following Hartlap et al. (2007), 
we then correct for bias in the inverse covariance matrix by 
multiplying it by the factor (A' im — Abm — iV-W'im- Below, 
we quote the significance level in terms of the fraction of 
simulations with a larger y 2 value than the observed map. 


C N {0 1 ,..., 0 2N -3) = (T(ni) • • • T{n N )), 
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Table 4. Probabilities of obtaining values for the \ 2 statistic 
of the A-point functions for the Planck fiducial ACDM model 
at least as large as the observed values of the statistic for the 
Planck 2015 temperature CMB maps with resolution parameter 
iV S ide = 64, estimated using the Commander, NILC, SEVEM, and 
SMICA methods. 


Function 

Comm, 

Probability [%] 

. NILC SEVEM SMICA 

2-pt. 

, . . . 97.2 

98.9 

97.4 

98.1 

Pseudo-coil. 3-pt. 

_ 92.1 

94.7 

91.8 

92.2 

Equil. 3-pt. 

, . . . 74.0 

80.4 

75.8 

79.0 

Rhombic 4-pt. 

, . . . 64.6 

70.9 

65.6 

65.9 


We analyse the CMB estimates at a resolution of 
7V S ide = 64, this being constrained by computational limita¬ 
tions. The results are presented in Fig. 2, where we compare 
the A-point functions for the data and the mean values es¬ 
timated from 1000 MC simulations. The probabilities of 
obtaining values of the x 2 statistic for the Planck fiducial 
ACDM model at least as large as the observed values are 
given in Table 4. 

It is worth noting that the values of the A-point func¬ 
tions for different angular separations are strongly corre¬ 
lated, and for this reason the figures show only one profile 
of each function in multi-dimensional space. Since the esti¬ 
mated probabilities take into account the correlations, they 
provide more reliable information on the goodness of fit be¬ 
tween the data and a given model than simple inspection 
of the figures. 

The results show excellent consistency between the 
CMB maps estimated using the different component- 
separation methods. No statistically significant deviations 
of the CMB maps from Gaussianity are found. Indeed, the 
slight preference for super-Gaussianity of the equilateral 3- 
point and 4-point functions observed for the 2013 data is 
now less pronounced. That may be caused by differences 
between the masks used for the analysis. Interestingly, the 
2-point function shows clear evidence of a lack of struc¬ 
ture for large separation angles. Such behaviour was origi¬ 
nally noted for the WMAP first-year data by Bennett et al. 
(2003), and has subsequently been discussed at length in 
the literature (Efstathiou 2004; Copi et al. 2007; Copi et al. 
2013). We will return to this issue in Sect. 5.2. 

4.4. Minkowski functionals 

Table 5. Probability P (x 2 > Xpianck) as a function of resolu¬ 
tion for the unnormalized, classical Minkowski functionals. 


Probability [%\ 


iV si de Comm. NILC SEVEM SMICA 

1024 . 91.4 90.7 95.5 95.8 

512 . 95.4 90.9 62.6 92.6 

256 . 55.8 34.5 55.9 55.9 

128 . 43.6 56.4 19.9 19.2 

64 . 59.3 37.8 22.7 80.0 

32 . 62.0 16.2 29.9 67.0 

16 . 43.4 45.8 47.7 31.0 


The Minkowski functionals (hereafter MFs) describe the 
morphology of fields in any dimension and have long been 


Table 6. Probability P (x 2 > Xpianck) as a function of resolu¬ 
tion determined using normalized MFs. 


Probability [%] 


A si de Comm. NILC SEVEM SMICA 

2048 . 97.2 77.7 99.0 93.0 

1024 . 93.1 98.0 90.2 92.6 

512 . 53.7 36.7 30.4 77.6 

256 . 89.0 85.9 96.8 58.1 

128 . 93.0 63.5 94.1 37.1 

64 . 37.1 70.4 54.1 62.5 

32 . 28.9 77.4 75.5 46.7 

16 . 33.1 39.4 44.1 38.8 


Table 7. Probability P (% 2 > Xpianck) as a function of needlet 
scale. 


Probability [%] 

Needlet scale {i range) Comm. NILC SEVEM SMICA 


3 (4,16). 32.1 36.1 40.4 39.8 

4 (8,32). 84.0 57.9 79.4 59.4 

5 (16,64) . 23.8 11.2 29.1 43.8 

6 (32,128) . 14.8 38.9 33.5 34.1 

7 (64,256) . 11.9 7.5 15.4 1.1 

8 (128,512) . 46.1 55.2 67.7 52.2 


used as estimators of non-Gaussianity and anisotropy in 
the CMB (see e.g., Mecke et al. 1994; Schmalzing & Buchert 
1997; Schmalzing & Gorski 1998; Komatsu et al. 2003; Erik- 
sen et al. 2004b; Curto et al. 2007; De Troia et al. 2007; 
Spergel et al. 2007; Curto et al. 2008; Hikage et al. 2008; 
Komatsu et al. 2009; Planck Collaboration XXIII 2014). 
They are additive for disjoint regions of the sky and invari¬ 
ant under rotations and translations. In the literature, the 
contours are traditionally defined by a threshold z/, usually 
given in units of the sky standard deviation (<7o). 

We compute MFs for the regions colder and hotter than 
a given threshold v. Thus, the three MFs, namely the area 
Vo( z/ ) = A(z/), the perimeter V\ (is) = C(z/), and the genus 
V 2 {y') = G(z/), are defined respectively as 

A 

V 0 (v)=A(v) = -^, (9) 

-^pix 

u (0 = c{v) = jl—XA ( 10 ) 

4A tot . 

V 2 (u) = G{v) = A— (JV hot - N cold ), (11) 

where N u is the number of pixels where AT /<tq > v, A p i X is 
the total number of available pixels, A tot is the total area of 
the available sky, Ah 0 t is the number of compact hot spots, 
A co id is the number of compact cold spots, and Si is the 
contour length of each hot spot. 

For a Gaussian random field in pixel space, the MFs can 
be written in terms of two functions: which depends 

only on the power spectrum, and Vk, which is a function 
only of the threshold v (see, e.g., Vanmarcke 1983; Pogosyan 
et al. 2009; Gay et al. 2012; Matsubara 2010; Fantaye et al. 
2015). The analytical expressions are 

14 (y) = A k Vk(v), (12) 
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Fig. 2. iV-point correlation functions determined from the N S id e = 64 Planck CMB 2015 temperature maps. Results are shown 
for the 2-point, pseudo-collapsed 3-point (upper left and right panels, respectively), equilateral 3-point, and connected rhombic 
4-point functions (lower left and right panels, respectively). The red dot-dot-dot-dashed, orange dashed, green dot-dashed, and 
blue long dashed lines correspond to the Commander, NILC, SEVEM, and SMICA maps, respectively. Note that the lines lie on top 
of each other. The black solid line indicates the mean determined from 1000 SMICA simulations. The shaded dark and light grey 
regions indicate the corresponding 68 % and 95 % confidence regions, respectively. See Sect. 4.3 for the definition of the separation 
angle 0. 


with 


derivative G\: 


Vkiy) = exp(-zv 2 /2)iJ fc _i(zv), k < 2, 

e - 2 

V3(iy) = erfc (v/ V2) ’ 


(13) A k 



1 L0 2 ( <7i \ k 

(27r)( fc +l)/ 2 U) 2 -k^k v \/2cro / ’ 

— ( ——L —) 2 

7T V \/2cr 0 / 


k < 2, 


(16) 

(17) 


and 

H n {v)=e'' 2 / 2 (-fAj l e-" 2 ' 2 . (15) 

The amplitude A^ depends only on the shape of the power 
spectrum Cn through the rms of the field ctq and its first 


where ujk = 7r k / 2 /T(k/2 + 1). 

Since this factorization is still valid in the weakly non- 
Gaussian case, we can use the normalized MFs, to focus 
on deviations from Gaussianity, with a reduced sensitivity 
to cosmic variance. 

Apart from the characterization of the MFs using full- 
resolution temperature sky maps, we also consider results 
at different angular scales. In this paper, two different ap¬ 
proaches are considered to study these degrees of freedom: 
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in real space via a standard Gaussian smoothing and degra¬ 
dation of the maps; and, for the first time, in harmonic 
space by using needlets. Such a complete investigation pro¬ 
vides an insight regarding the harmonic and spatial nature 
of possible non-Gaussian features detected with the MFs. 

First, we apply scale-dependent analyses in real space 
by considering the sky maps at different resolutions. The 
three classical MFs — area, contour length, and genus - 
are evaluated over the threshold range —3 < v < 3 in cr 
units, with a step of 0.5. This provides a total of 39 differ¬ 
ent statistics. The values of these statistics for the Planck 
data are all within the 95 % confidence region when com¬ 
pared with Gaussian simulations for all of the resolutions 
considered. A % 2 value is computed for each component- 
separation method by combining the 39 statistics and tak¬ 
ing into account their correlations (see e.g., Curto et al. 
2007, 2008). The corresponding covariance matrix is com¬ 
puted using 1000 simulations. The p- value of this y 2 test is 
presented in Table 5 for each component separation tech¬ 
nique and for map resolutions between 7V S id e = 1024 and 
Aside = 16. We find no significant deviations from Gaus- 
sianity for any of the resolutions considered. 

Then we consider the four normalized functionals de¬ 
scribed above. For every scale we used 26 thresholds ranging 
between —3.5 and 3.5 in a units, except for 6 = 640' where 
13 thresholds between —3.0 and 3.0 in a units were more 
appropriate. Table 6 indicates that no significant deviation 
from Gaussianity is found. 

Third, we tested MFs on needlet components. The 
needlet components of the CMB field as defined by Marin- 
ucci et al. (2008) and Baldi et al. (2009) are given by: 



Here, 77 (n) denotes the component at multipole £ of the 
CMB map T(n), i.e., 

T(n) = ^7Hn), (19) 

€ 

where ft G -S' 2 denotes the pointing direction, B is a fixed 
parameter (usually taken to be between 1 and 2) and 6(.) 
is a smooth function such that b 2 (£/Bi) = 1 for all £. 
Fantaye et al. (2015) show in a rigorous way that a general 
analytical expression for MFs at a given needlet scale j, 
which deals with an arbitrary mask and takes into account 
the spherical geometry of the sky, can be written as 

k 

VZ *Y, t V-i)Mvi, ( 20 ) 

i =0 


map and its first derivative given by 



Ct 


21+1 
47r ’ 


Ci 


21 + 1 i{t + 1 ) 

47T 2 


( 21 ) 

( 22 ) 


Implementing the MFs in needlet space has several ad¬ 
vantages: the needlet filter is localized in pixel space, hence 
the needlet component maps are minimally affected by 
masked regions, especially at high-frequency j; and the 
double-localization properties of needlets (in real and har¬ 
monic space) allow a much more precise, scale-by-scale, in- 
terpetation of any possible anomalies. While the behaviour 
of standard all-scale MFs is contaminated by the large cos¬ 
mic variance of the low multipoles, this is no longer the case 
for MFs evaluated at the highest needlet scales; in such cir¬ 
cumstances, the variance of normalized components may 
be shown to decrease steadily, entailing a much greater de¬ 
tection power in the presence of anomalies. Finally, and 
most importantly, the needlet MFs are more sensitive to 
the shape of the power spectrum than the corresponding 
all-scale MFs. 

The needlet parameters we use in this analysis are 
B = 2, j = 3,4, 5, 6, 7, 8. Since the masks in pixel space 
are map-resolution dependent, we also use different masks 
for each needlet scale. These new masks are constructed 
by multiplying the high-resolution common mask with the 
upgraded version of the appropriate low-resolution com¬ 
mon mask. For needlet scales j = 2 and j = 3, we use 
the common mask defined at 7V S ide = 16, and upgraded 
to A S ide = 2048. Similarly, for the higher needlet scales, 
j = 2 n , where n = 4, 5, 6, 7, 8, we use upgraded versions of 
the common masks defined at A S id e = 2 n . 

The results concerning needlet MFs from the 
Commander, NILC, SEVEM, and SMICA foreground-cleaned 
temperature maps for needlet scales B = 2, j = 4, 6,8 are 
shown in Fig. 3. All cases are computed using 26 thresholds 
ranging between —3.5 and 3.5 in a units. The figure shows 
the fractional difference between the Planck data and the 
FFP8 simulations in area (top panels), boundary length 
(middle panels), and genus (bottom panels) for different 
needlet scales. The j th needlet scale has compact support 
over the multipole ranges [2- 7-1 ,2 J+1 ]. All the scales 
we considered are consistent with the Gaussian FFP8 
simulations. This can be seen in Fig. 4, where we compare 
the data and simulation y 2 values, which are computed by 
combining the three MFs with an appropiate covariance 
matrix. The vertical lines in these figures represent the 
data, while the histogram shows the results for the 1000 
FFP8 simulations. We also show in Table 7 the p-values 
for the four component-separation methods, as well as 
all needlet scales we considered. Despite the relatively 
small p-values for some scales, the Planck temperature 
maps show no significant deviation from the Gaussian 
simulations up to £ max = 512, which corresponds to the 
maximum multipole of our highest-frequency needlet map. 


where to = 2, t\ = 0, and 1 2 = 47r are respectively the 
Euler-Poincare characteristic, boundary length, and area of 
the full sphere. The quantities Vk are the normalized MFs 
given in Eq. (13), while the needlet scale amplitudes A\ 
have a similar form as but with the variances of the 


4.5. Multiscale analyses 

Multiscale data analysis is a powerful approach for probing 
the fundamental hypotheses of the isotropy and Gaussian¬ 
ity of the CMB. The exploration of different scales (in an 
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Fig. 3. Needlet space MFs for Planck 2015 data using the four component-separated maps, Commander (red), NILC (orange), 
SEVEM (green), and SMICA (blue); the grey regions, from dark to light, correspond, respectively, to 1, 2, and 3cr confidence regions 
estimated from the 1000 FFP8 simulations processed by the Commander method. The columns from left to right correspond to the 
needlet parameters j = 4, 6, and 8, respectively; the j th needlet parameter has compact support over multipole ranges [2 J_1 , 2- 7+1 ]. 
The 4 = 2 J value indicates the central multipole of the corresponding needlet map. Note that to have the same range at all the 
needlet scales, the vertical axis has been multiplied by a factor that takes into account the steady decrease of the variance of the 
MFs as a function of scale. 


almost independent manner) not only helps to test the spe¬ 
cific predictions of a given scenario for the origin and evo¬ 
lution of the fluctuations, but also is an important check on 
the impact of systematic errors or other contaminants on 
the cosmological signal. 

There are several ways of performing a multiscale analy¬ 
sis, the simplest being to smooth/degrade the CMB map to 
different resolutions. However, in this section, we will focus 
on image processing techniques related to the application 
of wavelets and more general band-pass filtering kernels to 
the original CMB fluctuations. The advantage of wavelet¬ 
like analyses over scale degradation is clear: they allow the 
exploration of characteristics of the data that are related to 
specific angular scales. Wavelets have already been exten¬ 
sively used in the study of the Gaussianity and isotropy of 
the CMB (e.g., McEwen et al. 2007; Vielva 2007). Indeed, 
a wavelet-based (needlet) analysis of the Planck 2015 data 
has already been presented in Sect. 4.4. 

We recall that in the 2013 analysis, some of the applied 
estimators deviated from the null hypothesis. In particu¬ 
lar, it was determined that the cold area of the spherical 


Mexican hat wavelet (SMHW, Martinez-Gonzalez et al. 
2002) coefficients at scales of around 5° yielded a p-value 
of 0.3 %. In addition, we also found an excess in the kurto- 
sis of the wavelet coefficients on the same scales. Previous 
analyses (for a review, see Vielva 2010) have suggested that 
the “Cold Spot” (see Sect. 5.7) was the major contributor 
to these statistical outliers. 

In what follows, we will consider the application of the 
SMHW, together with matched filters for a 2D-Gaussian 
profile (GAUSS), and for generalized spherical Savitzky- 
Golay kernels (SSG, Savitzky & Golay 1964, see Ap¬ 
pendix A). 

The application of a filter !p(R,p) to a signal on the sky 
S(p) can be written as 

-bn ax m=£ 

L0 S ( R,P ) = E E semWf (R)Y em (p), (23) 

£=0 m=-l 

where p represents a given position/pixel, R parameterizes 
a characteristic scale for the filter (e.g., a wavelet scale), 
Wf ( R ) is the window function associated with the filter 
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Fig. 4. Histograms of % 2 for the Planck 2015 Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) foreground-cleaned 
maps analysed with the common mask. The % 2 is obtained by combining the three MFs in needlet space with an appropiate 
covariance matrix. The histograms are for the FFP8 simulations, while the vertical lines are for the data. The figures from left to 
right are for the needlet scales j = 4, 6, and 8, with the central multipoles £ c = shown in each panel. 


^(R,p), ^ max is the maximum multipole allowed by the cor¬ 
responding HEALPix pixelization, and Y^ m ( p ) is the spher¬ 
ical harmonic basis. Here, S£ m , the spherical harmonic co¬ 
efficients of the analysed map, are given by 

s im = JdflY; m (p)S(p), (24) 

where dU = d^sin^d^ and the asterisk denotes complex 
conjugation. Note that the filtered map (or the wavelet coef¬ 
ficient map, if ^{R,p) is a continous wavelet) conserves the 
statistical properties of the original map, since the convolu¬ 
tion is a linear operation. In particular, if S(p) is a Gaussian 
and statistically isotropic random signal, uis(R,p) is also 
Gaussian and statistically isotropic. 

In the present work, the signal S(p) corresponds to a 
temperature map T(p). Several statistics can then be com¬ 
puted from the derived filtered map as a function of the 
filter scale, in particular, the first moments (the dispersion 
or, the skewness Sr, and the kurtosis Kr ), the total area 
above/below a given threshold, and the peak distribution. 
These statistics are compared to the corresponding results 
determined from the FFP8 simulations to establish the de¬ 
gree of compatibility with the null hypothesis. 

4.5.1. First-order moments of the multiscale maps 

For the three filters considered (SMHW, GAUSS, and 
SSG84 3 ) the variance, skewness, and kurtosis are computed 
at 18 scales, R(arcmin) = {2, 4, 7, 14, 25, 50, 75, 100, 
150, 200, 250, 300, 400, 500, 600, 750, 900, 1050}. These 
scales are chosen to be consistent with previous analyses. 
They cover a wide angular range, and are selected so that 
the intervals between them increase with scale. Notice that, 
for a given scale, the three filters do not cover exactly the 
same multipole range, since that depends on the specific 
filter definition. This can be seen in Fig. 5: the SMHW is 
the narrowest filter, followed by SSG84, then GAUSS. The 
three filters have an equivalent effective £ max , but differ in 
the effective Anin- Overall, the differences between the fil¬ 
ters become smaller with increasing effective scale. In this 
paper, we refer to both the scale, R, and FWHM as pa¬ 
rameters defining the size of the filters. For the SMHW, 

3 The digits 8 and 4 denote the order of the spherical Savitzky- 
Golay kernel and the smoothing weight, described in Ap¬ 
pendix A. 


these are related by FWHM = R \/8 In 2, whereas for the 
GAUSS and SSG84 filters, the scale is defined to be half 
the FWHM. The latter definition is appropriate for filters 
that include pre-whitening since it is simple yet matches 
the Gspace bandwidth reasonably well. 

Following the procedure explained in PCIS13, after con¬ 
volution with a given filter, the common mask is extended 
to omit pixels from the analysis that could be contaminated 
by the mask. These pixels introduce an extra correlation 
between the data and the simulations, degrading the sta¬ 
tistical power of the comparison with the null hypothesis 
(see, e.g., Vielva et al. 2004). For a given scale R, the ex¬ 
clusion mask is defined by extending an auxiliary mask by 
a distance 2 R from its border, where the auxiliary mask is 
that part of the common mask related to residual diffuse 
Galacic emission (i.e., the auxiliary mask does not mask 
point sources). 

The following figures represent the upper-tail probabil¬ 
ity (UTP) for a given statistic, i.e., the fraction of sim¬ 
ulations that yield a value equal to or greater than that 
obtained for the data. In fact, as explained in PCIS13, if a 
given UTP is larger than 0.5, a new quantity is defined as 
mUTP = 1 — UTP. Therefore, mUTP is constrained to lie 
between 1/N and 0.5, where N is the number of simulations 
used for each statistic. 

Fig. 6 presents the comparison of the CMB temperature 
maps with the corresponding simulations for the SMHW, 
GAUSS, and SSG84 filters. The full mission Planck data 
confirm the results already obtained with the 2013 release 
for temperature. In particular, for the SMHW, we find 

(i) an excess of kurtosis (« 0.8%) at scales of around 300'; 

(ii) that the dispersion of the wavelet coefficients at these 
scales and at around 700' is relatively low (~ 1%); and 

(iii) that the dispersion of the wavelet coefficients at scales 
below 5' is significantly high (< 0.1 %). 

The excess of kurtosis has been previously associated 
with the “Cold Spot” (e.g., Vielva et al. 2004), and the 
low value of the standard deviation of the coefficients on 
large scales could be related to the low variance discussed 
in Sect. 5.1. Regarding the large dispersion of the coeffi¬ 
cients on the smallest scales, this can be understood either 
by the presence of residual foreground contributions (extra- 
galactic point sources) or by incomplete characterization of 
the true instrumental noise properties by the FFP8 simu¬ 
lations. We explore these possibilities with two additional 
tests undertaken with the SMHW. 
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Fig. 5. Comparison of the window functions (normalized to 
have equal area) for the SMHW (blue), GAUSS (yellow), and 
SSG84 (magenta) filters. The scales shown are 25 x (top) and 250 x 
(bottom). 


Figure 7 shows the significance of the statistics derived 
from the SEVEM-100, SEVEM-143, and SEVEM-217 maps. 
The three cleaned maps yield very consistent values of the 
mUTP for the standard deviation, skewness, and kurto- 
sis of the wavelet coefficients, with only small differences 
seen at small scales. This frequency-independence of the 
results argues against the foreground residuals hypothesis. 
Figure 8 presents the same statistics as applied to an es¬ 
timator of the noise properties of the CMB maps. This is 
derived from the half-difference of the half-ring data sets, 
which provides the best estimate of the noise properties 
of the full mission data set. However, since there is still a 
known mismatch in noise properties, any conclusions will be 
more qualitative than quantitative. Nevertheless, the noise 
study reveals that, at the smallest scales, there are some 
discrepancies with the FFP8 simulations, and in particular 
the estimated dispersion of the SMHW noise coefficients is 
higher than predicted. 


4.5.2. The area above/below a threshold 

In the context of the study of the Cold Spot, the area 
above/below a given threshold, as a function of the SMHW 
wavelet scale, has been demonstrated to provide a useful 
and robust statistic (e.g., Cruz et al. 2005), since it is rather 
independent of any masking required. Our previous analy¬ 
sis (PCIS13) confirmed that the CMB temperature fluctu¬ 
ations exhibit an anomalously large cold area on scales of 
around 10°, which can be mostly associated with the Cold 
Spot. Here, we extend the analysis by including results de¬ 
rived using the GAUSS and SSG84 filters. 

At a given scale R and threshold z/, the cold ( A and 
hot {Aft") areas of a filtered map are defined as 

a ~r = #{us(R,p) < -Z'}, (25) 

A+" = #{us(R,p)<+v}, (26) 

where the operator # represents the number of pixels p in 
which the condition defined between the braces is satisfied. 

Table 8 summarizes the results for the hot and cold ar¬ 
eas determined for the four CMB temperature maps anal¬ 
ysed with the common mask (and its associated exclu¬ 
sion masks). The results are similar to those obtained in 
2013, with some small differences on those scales related 
to the Cold Spot (between 200' and 400'). Specifically, the 
cold area is slightly less significant for smaller values of i?, 
whereas the anomalous behaviour remains for larger filter 
scales. The three filters are in reasonable agreement, but, as 
expected from Fig. 6, the SMHW yields higher significance 
levels than the SSG84 and GAUSS filters. However, it is 
worth recalling that, for a given scale, the three filters are 
not probing exactly the same multipole range and therefore 
some differences should be expected. 

In Fig. 9 we plot the areas for thresholds v > 3.0, where 
the threshold is defined in units of or, as determined from 
the SEVEM temperature map. The results for Commander, 
NILC, and SMICA are in good agreement with these. The 
panels refer to SMHW scales of R = 200', 250', 300', and 
400'. The most extreme value (in terms of ctr) for each area 
is indicated. 

The coldest area corresponds to the Cold Spot with the 
minimum value of the wavelet coefficient at the position 
(209°, —57°) in Galactic coordinates. The hottest area has 
already been identified in the WMAP data (e.g., Vielva 
et al. 2007). The results are insensitive to the choice of 
CMB temperature map that is adopted. It is clear that the 
southern Galactic hemisphere yields more anomalous sig¬ 
natures than the northern one. These results confirm the 
importance of the Cold Spot as the most extreme feature 
in the analysed sky. More insights about its nature are pro¬ 
vided in Sect. 5.7. 

4.5.3. Peak statistics 

The statistical properties of local extrema (both minima 
and maxima, which we refer to collectively as “peaks”) 
provide an alternative approach to search for evidence of 
non-Gaussianity in the data. Such peaks, defined as pix¬ 
els whose amplitudes are higher or lower than the corre¬ 
sponding values for all of their nearest neighbours, trace 
topological properties of the data. Peak locations and am¬ 
plitudes, and various derived quantities, such as their corre¬ 
lation functions, have previously been used to characterize 
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Fig. 6. Modified upper tail probabilities (mUTP) obtained from the analyses of the filter coefficients as a function of the filter scale 
R for the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) sky maps. From left to right, the panels correspond 
to standard deviation, skewness, and kurtosis results, when determined using the SMHW (top), GAUSS (middle), and SSG84 
(bottom) filters. The squares represent UTP values above 0.5, whereas circles represent UTP values below 0.5. 



Fig. 7. Modified upper tail probabilities (mUTP) obtained from the analyses of the SMHW coefficients as a function of the wavelet 
scale R for the SEVEM-100 (blue), SEVEM-143 (yellow), SEVEM-217 (magenta), and SEVEM (green) maps. From left to right, the panels 
correspond to the standard deviation, skewness, and kurtosis. 


the WMAP sky maps by Larson & Wandelt (2004, 2005) The statistical properties of peaks for a statistically 
and Hou et al. (2009). isotropic Gaussian random field were derived by Bond & 
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Fig. 8. Modified upper tail probabilities (mUTP) obtained from the analyses of the SMHW coefficients as a function of the wavelet 
scale R for the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) half-ring half-difference noise estimates. From left 
to right, the panels correspond to the standard deviation, skewness, and kurtosis. 


Table 8. Modified upper tail probability (mUTP ) for the cold 
(top) and hot (bottom) areas. Results are given for the is > 
4 ctr threshold of the SMHW, GAUSS, and SSG84 coefficients. 
The four most significant scales related to the Cold Spot feature 
are shown. An ellipsis (...) indicates that no area above that 
threshold was found in the data. 



Probability [%] 

Area 

Scale Comm. NILC SEVEM SMICA 
[arcmin] 



SMHW 





200 3.8 

5.1 

3.7 

3.8 

Cold .... 

250 1.4 

2.4 

1.4 

1.4 


300 0.4 

1.5 

0.4 

0.4 


400 0.9 

0.9 

0.9 

0.9 


200 2.0 

2.6 

1.7 

1.5 

Hot .... 

250 2.4 

3.0 

2.1 

2.0 


300 4.2 

400 

5.0 

4.1 

3.9 


GAUSS 





200 1.7 

2.7 

1.7 

1.7 

Cold .... 

250 1.2 

1.2 

1.2 

1.2 


300 1.6 

400 

1.8 

1.2 

1.8 


200 2.9 

3.5 

2.8 

2.6 

Hot .... 

250 5.7 

300 

400 

6.4 

5.6 

5.4 


SSG84 





200 9.4 

11.0 

9.4 

9.0 

Cold .... 

250 12.3 

13.4 

10.8 

12.3 


300 1.4 

2.6 

1.4 

1.5 


400 0.9 

1.9 

0.8 

0.9 


200 1.1 

1.8 

1.0 

0.9 

Hot .... 

250 4.8 

300 

400 

5.1 

4.5 

4.3 


Efstathiou (1987). The integrated number density of peaks, 
7i p k (composed of maxima and minima with corresponding 
densities n max and n m i n ), with amplitudes x above a certain 


threshold is = x/g is given by 



where cr is the rms fluctuation amplitude measured on the 
sky, and 7 is the spectral shape parameter of the underly¬ 
ing field. Uncharacteristically cold and hot spots are then 
manifested as extreme outliers in the peak values, and can 
constitute evidence for non-Gaussianity or deviation from 
isotropy. 

Here, we consider the peak statistics of the Planck 
component-separated temperature maps at 7V S id e = 2048. 
The maps are pre-whitened as described in Appendix A. 
This step allows the construction of an estimator that is 
nearly optimal with respect to the fiducial CMB properties. 
After application of the common mask, weighted convolu¬ 
tions of the data are performed with either SSG or GAUSS 
kernels of variable scale. In order to avoid potential contam¬ 
ination by boundary effects, the mask is extended by reject¬ 
ing pixels with an effective convolution weight that differs 
from unity by more than 12 %. Peaks are extracted from 
the filtered map (removing any that are adjacent to masked 
pixels), their positions and values are recorded for further 
analysis, and their cumulative density function (CDF) is 
constructed by sorting peak values. Table 9 presents peak 
counts for the component-separated sky maps for several 
different kernels and representative filtering scales, together 
with the number of peaks that are common to all maps. 
There is excellent agreement between the various CMB es¬ 
timates. All statistical inference is then performed by com¬ 
parison of the peak distributions derived from the data with 
equivalently processed simulations. As an internal consis¬ 
tency check, the properties of the FFP8 simulations are 
found to be in agreement with the predictions of Eq. (27). 

Figure 10 presents the distributions of peaks for the 
SMICA CMB map filtered with two representative kernels 
on scales of 40' and 800' FWHM. The lower panels show 
the empirical peak CDFs as a function of peak value x, 
defined for a set of n peaks at values {X{\ as 

*.<*>-'**.-{«; • < 28 > 

i= 1 ^ 
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R = 200 arcmin 






Fig. 9. Cold and hot areas for thresholds v > 3.0 as determined 
from the SEVEM temperature map. From top to bottom, the maps 
are for SMHW scales of R = 200', R = 250', R = 300', and 
R = 400'. 


For plotting purposes alone, the horizontal axis is scaled in 
units of a defined by Eq. (27) and derived from the underly¬ 
ing median CDF, F(x), of the simulations. The upper pan¬ 
els show the difference between the observed and median 
simulated CDF values, yjn[F n (x ) — F(x)], with the grey 
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Fig. 10. Cumulative density function of the peak distribu¬ 
tion for the SMICA CMB temperature map. The top row shows 
the peak CDF for maps filtered with a GAUSS kernel of 40 / 
FWHM. The bottom row shows the corresponding peak CDF 
for an SSG84 kernel of 800 7 FWHM. The spectral shape parame¬ 
ter 7 (see Eq. 27) is the best-fit value for the simulated ensemble, 
as indicated by the cyan circle in Fig. 11. Similar results are ob¬ 
tained for the other component-separation methods. 


bands representing the 68.3%, 95.4%, and 99.7% regions 
of the simulated CDF distributions. The maximal value of 
this difference defines a Kolmogorov-Smirnov (KS) devia¬ 
tion estimator: 

K n = \fn sup | F n (x) - F(x) | . (29) 

X 

This forms the basis of a standard KS test of consistency 
between the two distributions. Although the KS deviation 
has a known limiting distribution, we also derive its CDF 
directly from the simulations. 

The temperature peak distributions in Fig. 10 are con¬ 
sistent with Gaussian peak statistics, apart from a single 
anomalously cold peak on scales around 800' FWHM. This 
corresponds to the previously reported Cold Spot. Although 
this exercise confirms that the Cold Spot is a rare cold fea¬ 
ture, as already noted by Cruz et al. (2005) and confirmed 
in this paper, the most peculiar characteristic of the Cold 
Spot is not its coldness, but rather its size. A more detailed 
analysis of its nature is presented in Sect. 5.7. 
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a 



Fig. 11. Distribution of best-fit Gaussian peak CDF spectral 
shape parameters, a and 7 (as defined in Eq. 27), recovered 
from 1000 simulations, as indicated by the black dots and the 
smoothed density map and compared to those derived for the 
observed sky (shown by the red star). The blue contours enclose 
68 % and 95 % of the parameter distribution, and the cyan cir¬ 
cle represents the best-fit parameters for the median peak CDF 
determined from simulations. The upper panel shows the peak 
CDF parameters for the SMICA map filtered with a GAUSS ker¬ 
nel of 40 7 FWHM. The lower panel shows the corresponding 
peak CDF for an SSG84 kernel of 800' FWHM. Similar results 
are obtained for the other component-separation methods. 


The probability that the observed sky exceeds the value 
of the KS deviation for the adopted fiducial cosmology can 
be determined by counting the number of simulations with 
K n > > id sky) . The p -values for the KS test comparing the 
CDF of the observed sky with the median peak CDF de¬ 
rived from simulations for several different kernels and rep¬ 
resentative scales are summarized in Table 10. The similarly 
derived p -values for the total peak counts are summarized 
in Table 11. Most of the results indicate that the two distri¬ 
butions are highly consistent, with the exception of results 
for the SSG84 filter on scales of about 500' FWHM. This 
deviation appears to be related to a hemispherical asym¬ 
metry in the peak CDFs, and will be discussed further in 
Sect. 5.6. 



Mean amplitude of extrema 
o 



Fig. 12. Cumulative density function of the mean amplitude 
of all extrema, maxima (red) and minima (blue), derived from 
simulations, compared to the equivalent values observed for the 
SMICA CMB temperature map. The upper panel shows the peak 
mean amplitudes for maps filtered with a GAUSS kernel of 4(fi 
FWHM. The lower panel shows the corresponding peak CDF 
for an SSG84 kernel of 800 x FWHM. Similar results are obtained 
for the other component separation methods. Since the filter 
kernel normalization is free, and the pre-whitened map to which 
the filter is applied is dimensionless, the plots are essentially in 
arbitrary units. 

Table 9. Peak counts in maps filtered to different scales. 


Number of minima/maxima 

Filter Scale Comm. NILC SEVEM SMICA Match 
[arcmin] 


SMHW 

200. 176/187 170/178 173/182 169/182 161/169 

250 . 105/105 104/103 107/123 105/107 97/ 99 

300 . 70/ 70 71/ 70 70/ 72 68 / 71 66 / 66 

400 . 43/ 32 46/ 32 44/ 31 43/ 33 37/ 30 

GAUSS 

200 . 152/170 152/166 157/179 156/165 142/155 

250 . 102/ 93 104/ 95 108/ 99 99/101 92/ 85 

300 . 60/ 63 57/ 62 63/ 64 56/ 62 50/ 53 

400 . 33/ 28 29/ 29 31/ 33 29/ 28 24/ 27 

SSG84 

200 . 180/187 178/183 180/185 183/183 167/175 

250. 131/119 118/114 122/123 121/110 109/103 

300 . 68 / 69 73/ 68 73/ 73 70/ 68 56/ 61 

400 . 29/ 35 29/ 36 29/ 32 30/ 38 27/ 27 


One can also test whether the observed values of the 
parameters, a and 7 as defined in Eq. (27), are consistent 
with the simulation ensemble, under the assumption that 
the peak distributions in the Planck data are described by 
a Gaussian peak CDF. Figure 11 demonstrates the consis¬ 
tency of the best-fit values of these parameters, correspond¬ 
ing to the peak distributions in Fig. 10, with equivalent 
values derived from the simulations. 
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Fig. 13. Fraction of the Gaussian random field realizations in which the coldest peak is as cold as or colder than that observed, 
as a function of SMHW filter scale for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). 


Table 10. Modified upper tail probability (mUTP) for the KS 
test, comparing the data with the median peak CDF derived 
from simulations. 


Filter Sc 
[arcmii 

Probability [%] 

zale Comm. 

i] 

NILC 

SEVEM 

SMICA 


SMHW 



200 . 

_ 22.0 

42.8 

45.9 

40.5 

250 . 

.... 11.3 

17.6 

3.1 

11.4 

300 . 

_ 49.4 

38.5 

38.4 

32.1 

400 . 

.... 32.6 

24.7 

35.3 

24.7 


GAUSS 



200 . 

_ 41.3 

46.6 

14.4 

47.2 

250 . 

_ 43.7 

34.8 

7.6 

48.4 

300 . 

.... 46.3 

9.9 

28.0 

7.7 

400 . 

_ 30.7 

5.6 

35.8 

6.6 


SSG84 



200 . 

.... 37.1 

36.7 

24.0 

37.5 

250 . 

_ 0.5 

1.7 

0.8 

5.4 

300 . 

_ 17.5 

12.2 

0.3 

9.3 

400 . 

47.4 

44.6 

47.5 

47.8 


Inspired by the analysis of the WMAP first-year data 
in Larson & Wandelt (2004) which found fewer extreme 
peaks than expected, we additionally evaluate whether the 
distributions of maxima and minima are separately con¬ 
sistent with simulations. The mean of all maxima, and 
the negative of the mean of all minima, are calculated for 
the filtered map, and the observed values are compared to 
the simulated distributions in Fig. 12. The observed min¬ 
ima/maxima means are found to be in good agreement with 
the fiducial values. 

The probability that the coldest peak seen on the sky is 
consistent with the adopted fiducial cosmology is evaluated 
as a function of both filter shape and size by counting the 

number of simulations with boldest < boldest- The resu ^ s 
obtained for the SMHW filter are summarized in Fig. 13. 
Consistent behaviour is seen when the GAUSS and SSG84 
filters are applied. The error bars represent the sampling 
uncertainty due to the finite number of realizations, and 


Table 11. Modified upper tail probability (mUTP) for the to¬ 
tal peak count, comparing the data with the peak count CDF 
derived from simulations. 


Filter Sc 
[arcmii 

Probability [%] 

ale Comm. 

i] 

NILC 

SEVEM 

SMICA 


SMHW 



200 . 

_ 6.1 

36.9 

16.2 

27.2 

250 . 

.... 32.9 

47.5 

1.0 

25.6 

300 . 

_ 48.8 

51.7 

44.7 

44.3 

400 . 

.... 33.8 

16.2 

34.6 

26.4 


GAUSS 



200 . 

7.1 

11.2 

0.7 

8.7 

250 . 

_ 18.2 

11.2 

2.1 

8.1 

300 . 

.... 29.0 

12.8 

48.2 

10.0 

400 . 

_ 11.9 

3.0 

26.6 

2.8 


SSG84 



200 . 

0.2 

3.0 

1.0 

1.7 

250 . 

_ 0.1 

1.7 

0.1 

2.1 

300 . 

9.3 

22.6 

50.7 

12.2 

400 . 

0.3 

0.4 

0.1 

2.3 


are determined using a bootstrap method. As the filters 
overlap substantially, different points are highly correlated. 
The Planck CMB maps are consistent with the expecta¬ 
tions of a statistically isotropic Gaussian model. The most 
significant deviation, found at an effective filter bandwidth 
given by i = 20, is attributable to a single region on the 
sky — the Cold Spot. 

4.5.4. Peak locations as a function of scale 

The application of a filter kernel of variable size to a map 
extends it into what can be considered a “multiscale space,” 
such that features on different scales are represented by a 
one-parameter family of smoothed maps. This technique 
is often used for feature detection and mathematical mor¬ 
phology analysis. Here, we introduce a morphological de¬ 
scription of temperature maps based on the peak connect¬ 
edness graph in multiscale space, and apply this technique 
to a statistical analysis of the Planck CMB data. Like most 
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Fig. 14. Peak positions and CDF rank summarized for all filtering scales. The three sky-view panels in the top row show a 
Lambert projection of the north pole, the usual full sky Mollweide projection, and a Lambert projection of the south pole. The 
lower panel shows the peak heights (in percentile of the peak distribution on the horizontal axis) as a function of filter scale (on 
the vertical axis, in logarithmic scale), truncated to larger scales for clarity. Circles represent peaks (nodes of the graph) coloured 
according to their percentile level, and scaled according to kernel size. Black lines represent edges connecting peaks at different 
scales (according to a minimal distance measure). The components connected to the coldest and hottest peaks at any scale are 
highlighted by thicker edges, and are navy blue and dark red in colour. Note that there are thick lines that do not touch the 0 and 
1 percentiles in the plot view. Those edges are connected to extreme percentile values, but at scales smaller than those shown in 
the plot. The Cold Spot is represented by the connected nodes that have the smallest percentiles except for the coarsest scale in 
the plot view. 


morphological analyses, it is equally applicable to intrinsi¬ 
cally non-Gaussian maps, but here we focus on the Gaussian 
random field statistics and attempt to understand what fea¬ 
tures of the CMB temperature map are responsible for the 
Cold Spot. 

To construct a multiscale representation, we trace the 
location of the peaks in the smoothed, whitened CMB 
map as the smoothing scale is varied. As the smoothing 
scale increases, peaks merge and the total peak count de¬ 
creases. Linking closest peak neighbours in position-scale 
space, from the finest to the coarsest resolution, produces 
an acyclic graph that encapsulates the peak “merger tree” 
history as the scale is varied. A summary of all the peak po¬ 
sitions and CDF ranks for the SSG84 filter kernel on scales 
ranging from 120' to 1200' FWHM is shown in Fig. 14. The 
peaks are represented by discs of varying size (reflecting 
the filter scale) and colour (reflecting the peak tempera¬ 
ture rank), with peaks at all scales projected onto a single 
map. The lower panel shows the peak linkage graph on the 
coarser scales; for the statistical analysis 81 filter scales are 
used, log-spaced from 120' to 1200'. Peaks of the same type 
(i.e., maxima to maxima and minima to minima) are linked 
to the closest peak on the coarser scale according to a dis¬ 
tance measure, d s 2 + d/ 2 , where ds 2 is the metric on the 
unit sphere, and d/ 2 is the difference of peak temperature 



Fig. 15. Distribution of node degrees in the multiscale peak 
linkage graph determined for the SMICA map (cyan), compared 
with the median (red line), first to third quartile (blue box), and 
95 % (whiskers) derived from 1000 FFP8 simulations. 


ranks (but only if that distance is within a predetermined 
fraction of the filter scale FWHM). 

The resulting peak linkage graph is then analysed for 
connectedness. The simplest quantifiable measure is the 
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node-degree distribution, shown in Fig. 15 for SMICA. The 
node-degree distribution is highly peaked at 2; this pop¬ 
ulation corresponds to a single peak being traced across 
multiple scales. Pre-whitening effectively decorrelates the 
Gaussian map across different scales, so that the resulting 
node distribution shows a sizeable population of degree 0 
and 1 nodes. When compared to the linkage graphs derived 
from the simulation ensemble, the node-degree distribution 
of the peak linkage graph derived from Planck CMB data 
is consistent, with a slight excess in node counts of degrees 
5 and 6. 

5. Anomalies in the microwave sky 

The previous section established the lack of evidence for 
significant non-Gaussianity in the Planck data. Here we 
consider several important anomalies that were originally 
detected in the WMAP sky maps, and later confirmed in 
the analyses described in PCIS13. Many of these are con¬ 
nected to evidence for a violation of isotropy, or to a pre¬ 
ferred direction, in the CMB. Tests that involve dipolar 
power asymmetry, either directly or via measures of direc¬ 
tionality, are collected together in Sect. 6. In this section 
we consider only those anomalies not directly related to 
dipolar power asymmetry. 

The microwave sky is intrinsically statistically 
anisotropic due to our motion with respect to the 
CMB rest frame. The resulting Doppler boosting effect, 
introduced in Sect. 1, was detected in the 2013 Planck data 
(Planck Collaboration XXVII 2014). For completeness, 
Appendix B repeats the analysis with the Planck full 
mission data set. However, since the effects of Doppler 
boosting are now included in the simulations used for 
that analysis, this constitutes a consistency check for 
this release. More importantly, since both the data and 
simulations now include the effect, it is not necessary 
to consider deboosted data in many of the studies re¬ 
ported here, unlike in PCIS13 (although one exception in 
Sect. 6.4 makes use of unboosted simulations to search 
for the frequency-dependent signature of the effect in the 
SEVEM-100, SEVEM-143, and SEVEM-217 sky maps). 

Before presenting our results, we return to the issue of 
a posteriori correction, which particle physicists refer to 
as correcting for the “look-elsewhere effect” (LEE). Since 
there are many tests that can be performed on the data to 
look for a violation of statistical isotropy, we expect some 
to indicate detections at, for example, roughly 3 cr levels, 
since even a statistically isotropic CMB sky is a realization 
of an underlying statistical process corresponding to many 
independent random variables. However, in the absence of 
an existing theoretical framework (i.e., a physical model) 
to predict such anomalies, it is difficult to interpret their 
significance. It is then necessary, and equally challenging, 
to address the question of how often such detections would 
be found for statistically isotropic Gaussian skies. Unfortu¬ 
nately, it is not always clear how to answer this question. 

There will always be a degree of subjectivity when de¬ 
ciding exactly how to assess the significance of these types 
of features in the data. As an example, one could argue 
that the large-scale dipole modulation signal we see is com¬ 
ing specifically from super-Hubble modes, in which case 
performing a LEE correction for dipole modulation that 
could have been seen on small scales (£ > 100) would not 
make sense. Models for such a super-Hubble modulation ex¬ 


ist and an example was examined in Planck Collaboration 
XX (2015), the conclusion being that the model could only 
explain part of the dipole modulation and that the allowed 
part was perfectly consistent with cosmic variance. 

In this paper, we adopt a pragmatic approach. When 
there is a clear mechanism for doing so, we attempt to cor¬ 
rect for the “multiplicity of tests,” or the possible ways in 
which an anomalous signal might have been detected but 
was not, as a consequence of any a posteriori (data-driven) 
choices made in searching for it. In such cases, a strong 
dependence of the significance on the correction would in¬ 
dicate that we should be cautious about the uncorrected 
result. When such an obvious correction is not possible, we 
clearly describe the methodology applied to the data and 
its limitations. With this approach, we also recognize that 
any statistical assessment is partially subjective, including 
those that purport to correct for the LEE. 

Although many of the observed effects described in this 
and the next section may elude theoretical prediction today, 
we continue to highlight them since there is a real possibil¬ 
ity that the significance of one or more might increase at a 
later date, perhaps when polarization data are included in 
the analysis, and lead to new insights into early Universe 
physics. Alternatively, such observations may directly moti¬ 
vate the construction of models that can make predictions 
for features that can be sought in new data sets. This is 
particularly the case for anomalies on the largest angular 
scales, which may have a specific connection to inflation. 

5.1. Variance , skewness , kurtosis 

Previous analyses of the WMAP data (Monteserfn et al. 
2008; Cruz et al. 2011; Gruppuso et al. 2013) have reported 
that the variance of the CMB sky is lower than that of 
simulations based on the ACDM model. PCIS13 confirmed 
this, and proposed a possible explanation of the apparent 
incompatibility of the observed variance with a fiducial cos¬ 
mological model that has been determined from the same 
data set. Specifically, whilst the map-based variance is dom¬ 
inated by contributions from large angular scales on the sky, 
the cosmological parameter fits are relatively insensitive to 
these low-order C modes, and are instead largely dominated 
by scales corresponding to £ > 50. Therefore the variance 
of the map appears to be anomalous, since there is a dearth 
of large-angular-scale power compared to the model. 

In Sect. 4.1, we again confirmed the presence of low 
variance in the data. Here, we extend the analysis to inves¬ 
tigate which angular scales are responsible for the low vari¬ 
ance by applying the unit variance estimator to lower res¬ 
olution component-separated maps, specifically those from 
lV S ide = 1024 to A S ide = 16, with the corresponding com¬ 
mon masks, and then comparing the results with those de¬ 
termined from 1000 MC simulations. The results are shown 
in Fig. 16. 

All of the component-separation methods that we con¬ 
sider yield very consistent results which indicate an increas¬ 
ingly anomalous low variance at lower resolutions, with the 
lower-tail probability reaching a minimum value of 0.5 % 
at Aside = 10. We then consider the impact of a possible 
look-elsewhere effect by evaluating the minimum lower tail 
probability of each simulation irrespective of the A S id e reso¬ 
lution at which it occurs. By comparing the distribution of 
these values with that of the data, we infer that the prob¬ 
ability is slightly weakened to a value of about 1 %. These 
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Fig. 16. Lower tail probability of the variance (top panel), skew¬ 
ness (centre panel), and kurtosis (bottom panel) obtained at 
different resolutions from the Commander (red), NILC (orange), 
SEVEM (green), and SMICA (blue) sky maps. 


Table 12. Lower-tail probability for the variance, skewness, and 
kurtosis of the low resolution N S id e = 16 component-separated 
maps obtained with the common mask and two extended ver¬ 
sions thereof. 



Probability [%] 

Method 

Variance Skewness Kurtosis 


Common mask (/ s k y — 58 %) 


Commander . . 

0.5 

14.6 

88.4 

NILC. 

0.5 

16.9 

87.1 

SEVEM . 

0.5 

17.2 

84.8 

SMICA . 

0.5 

16.6 

82.7 


/ 3ky =48% 



Commander . . 

0.1 

29.4 

65.0 

NILC. 

0.1 

29.6 

60.8 

SEVEM . 

0.1 

29.4 

62.4 

SMICA . 

0.1 

29.4 

57.3 


/ sky = 40% 



Commander . 

0.4 

35.2 

32.4 

NILC . 

0.4 

34.4 

28.7 

SEVEM . 

0.4 

34.3 

30.2 

SMICA . 

0.4 

33.8 

25.5 


results are compatible with a lack of power on large angu¬ 
lar scales. Since the variance estimator is heavily weighted 
towards low-£ modes, this has an increasing impact on the 
observed variance when going from high to low resolution 
sky maps. Conversely, the skewness and kurtosis are consis¬ 
tent with the simulations, although there is some indication 
of a weak scale-dependence (albeit at low significance). 

We also investigate the stability of the results at N s ^ e = 
16 with respect to the possible presence of residual fore¬ 
grounds by considering two additional masks obtained by 
extending the edge of the Wide = 16 common mask by 5° 
and 9°, reducing the usable sky fraction from 58% to 48% 
and 40 %, respectively. We then re-apply the unit variance 
estimator to the low resolution component-separated maps 
with these masks and determine the variance, skewness, and 
kurtosis values (see Table 12). 


Table 13. Probabilities of obtaining values for the S ±/2 and 
Xo statistics for the Planck fiducial ACDM model at least as 
large as the observed values of the statistic for the Planck 2015 
temperature CMB maps with resolution parameter Wide = 64, 
estimated using the Commander, NILC, SEVEM, and SMICA maps. 
We show also the corresponding estimation of the global p -value 
for the S(x) statistic. 


Statistic 

Comm. 

Probability [%] 

NILC SEVEM 

SMICA 

Sl/ 2 . 

_ 99.5 

99.6 99.5 

99.6 

S(x) (global) . 

_ 97.7 

97.8 97.8 

97.9 

X 2 o(0 > 60 ) . 

_ 98.1 

98.8 98.1 

98.4 


The results from 48 % of the sky reveal that only 1 sim¬ 
ulation in 1000 is found to be more anomalous (i.e., exhibit 
lower variance) than the observed map. In addition, both 
the skewness and kurtosis become more compatible with the 
ACDM model. With the more aggressive mask, the lower- 
tail probability slightly increases again. However, given the 
limited number of pixels involved in the analysis, this shift 
may be related to the effects of sample variance. 

Overall, our results may be explained by the presence 
of a low-variance anomaly in the primordial CMB signal — 
the stability of the low-variance significance argues against 
foreground contamination being responsible for the lack 
of observed power. This is reinforced by the decrease in 
variance when regions close to the common mask borders, 
where foreground residuals are most likely to be observed, 
are omitted from the analysis. 


5.2. N-point correlation function anomalies 
5.2.1. Lack of large-angle correlations 

We first reassess the lack of correlation seen in the 2-point 
angular correlation function at large angular separations as 
reported in Sect. 4.3, and previously noted for both WMAP 
and the 2013 Planck temperature maps (Bennett et al. 
2003; Copi et al. 2013). We attempt to quantify this lack 
of structure using the statistic proposed by Spergel et al. 
(2003): 




d(cos$) , 


(30) 


where C^#) is our estimate of the 2-point correlation func¬ 
tion. Generally, the upper limit on the integral has been 
taken to correspond to a separation angle of 60°, possibly 
(as noted by Copi et al. 2009) motivated by the COBE- 
DMR 4-year results (Hinshaw et al. 1996). Inspection of 
the top panel of Fig. 2 suggests that the Planck 2-point 
function lies close to zero between 80° and 170°, but for 
consistency with previous work we compute the statistic 
Si/ 2 , for x = cos60° = The results are presented in Ta¬ 
ble 13. We find that the data indeed show a lack of correla¬ 
tions on large angular scales, with a significance consistent 
with that found by Copi et al. (2013) (although note that 
the sense of the p -values differs between the papers). 

Possible criticisms of the S i/ 2 statistic include that it 
has been designed a posteriori to test for a lack of large- 
angle correlations, and that it does not account for the 
high degree of correlation between bins at different angular 
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Table 14. Probabilities of obtaining values for the x 2 statistic 
for the Planck fiducial ACDM model at least as large as the 
observed values of the statistic for the Planck 2015 temperature 
CMB maps with resolution parameter A S id e = 64, estimated 
using the Commander, NILC, SEVEM, and SMICA maps. 



Probability [%] 


Statistic 

Comm. NILC SEVEM 

SMICA 

x 2 (d < 60°). 

. . . . 91.5 93.3 91.6 

91.7 

X 2 (0 > 60°) . 

_ 96.8 98.3 96.9 

98.1 


scales. We can address these concerns, at least in part, by 
considering a modified version of the commonly used and 
well understood x 2 statistic used in previous studies. In or¬ 
der to test the same hypothesis as the Si/ 2 statistic — that 
there are no correlations on scales larger than some angular 
cut-off — we do not subtract an averaged 2-point function 
when computing the x 2 > i.e., we use a statistic defined as 

xl(o min ,e max )= 2 c 2 {e i )u- l c 2 {e J ) , (3i) 

b 3 ~*min 

where i m i n , i max denote the index of the bins corresponding 
to the minimum and maximum value of the separation an¬ 
gles 6 m i n and # max , respectively. In this analysis, we adopt 
#min = 60° and # max = 180°. M^- is the covariance ma¬ 
trix given by Eq. (8), estimated using MC simulations cor¬ 
responding to the fiducial ACDM model. The results are 
shown in Table 13. The significance level of the anomaly is 
slightly smaller for the Xo statistic than that derived with 
*S ]_/2 • We note that this statistic is closely related to the 
A{x) measure proposed by Hajian (2007). 

A further potential criticism of the S' 1 / 2 statistic relates 
to the a posteriori choice of 60° for the separation angle 
that delineates the interesting region of behaviour of the 
correlation function. We therefore consider the generalized 
statistic S(x) and compute its value for all values of x, both 
for the data and for the simulations. Then, for each value 
of x, we determine the number of simulations with a higher 
value of S(x), and hence infer the most significant value of 
the statistic and the separation angle that it corresponds 
to. However, since such an analysis is sensitive to the LEE, 
we define a global statistic to evaluate the true significance 
of the result. Specifically, we repeat the procedure for each 
simulation, and search for the largest probability irrespec¬ 
tive of the value of x at which it occurs. The fraction of 
these probabilities higher than the maximum probability 
found for the data defines a global p-value. As seen in Ta¬ 
ble 13, this corresponds to values of order 98% for all of 
the CMB estimates. 

The previous analyses essentially test how consistent the 
observed 2-point correlation function data is with a lack of 
correlations on large angular scales, in particular for separa¬ 
tion angles 0 > 60°. A conventional x 2 statistic allows us to 
test the consistency of this quantity with the predictions of 
the ACDM model. In this case, the statistic is defined as in 
Eq. (7), except that we constrain the computations to those 
bins that correspond to the intervals defined by 6 < 60° and 
0 > 60°. The results of these studies are shown in Table 14. 

The analysis for 0 < 60° indicates that the observed 2- 
point function is a good match to the mean 2-point function 
predicted by the ACDM model. Moreover, for 0 > 60° the 


results suggest that the problem is that the fit of the data 
to the model is too good, and this is even more pronounced 
for an analysis in the full separation angle range. 

Overall, the tests indicate an unusually good fit of the 
observed 2-point function both to zero and to the predic¬ 
tions of the ACDM model for angles above 60°. This prob¬ 
lem may be related to the fact that the theoretical variance 
for the best-fit model is larger than the observed value at 
large scales, so that the simulations based on this model 
that have been used in all of the statistical tests may over¬ 
estimate the variance of the 2-point function. 

5.2.2. Hemispherical asymmetry 

We now turn to a reassessment of the asymmetry between 
the real-space A-point correlation functions computed on 
hemispheres reported previously for the WMAP (Eriksen 
et al. 2005) and Planck 2013 temperature maps (PCIS13). 
We initially focus the analysis on the hemispheres de¬ 
termined in the ecliptic coordinate frame for which the 
largest asymmetry was observed. However, we also carry 
out the corresponding calculations in other relevant refer¬ 
ence frames, such as those defined by the Doppler boost 
(DB, see Sect. 6.4, Appendix B, and Planck Collaboration 
XXVII 2014) and the dipole modulation (DM, see Sect. 6.2) 
directions. We use the same configurations of the A-point 
functions as described in Sect. 4.3. However, here the func¬ 
tions are not averaged over the full sky and depend on 
a choice of specific direction, so they constitute tools for 
studying statistical isotropy rather than non-Gaussianity 
(Ferreira & Magueijo 1997). 

As in Sect. 4.3, we analyse the CMB estimates at a res¬ 
olution of A^ide = 64 and quantify their agreement with the 
fiducial cosmological model using a x 2 statistic. The results 
determined from the Planck 2015 temperature data for the 
ecliptic hemispheres are shown in Fig. 17. If we consider 
that the x 2 statistic itself can act as a measure of fluc¬ 
tuation level, then asymmetry between the two measured 
hemispheres can be quantified by the ratio of the corre¬ 
sponding x 2 values. The probabilities of obtaining values 
of the x 2 statistic or ratio for the Planck fiducial ACDM 
model at least as large as the observed values are given in 
Table 15. Since we do not have any predictions concerning 
the behaviour of a given hemisphere, in the case of the x 2 
ratios we provide the complementary probabilities of the 
2-tailed statistic. 

The significance levels of the 3- and 4-point functions in 
the northern hemisphere are nominally very high, exceeding 
99.9 % for the pseudo-collapsed 3-point function. However, 
proper interpretation requires that one recognize that the 
analysis is affected by a posteriori choices for the smoothing 
scale and reference frame defining the hemispheres. This 
typically leads to an overestimation of the significance of 
the results. Accounting for such effects requires the repeti¬ 
tion of the analysis for all possible reference directions and 
also for data at other resolutions. Unfortunately, because 
of computational limitations, such an extended analysis is 
not possible for these higher-order statistics. Nevertheless, 
the observed properties of the Planck data are consistent 
with a clear lack of fluctuations in a direction towards the 
north ecliptic pole. However, the x 2 - r &tio statistic indicates 
a slightly smaller significance level for the asymmetry, not 
exceeding 99 % for any of the A-point functions. 
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Table 15. Probabilities of obtaining values for the x 2 statistic 
and ratio of x 2 of the N -point functions shown in Fig. 17 for the 
Planck fiducial ACDM model at least as large as the observed 
values of the statistic for the Planck 2015 CMB maps estimated 
on northern and southern ecliptic hemispheres. 


Probability [%] 


Hemisphere Comm. NILC SEVEM SMICA 

2-point function 

Northern. 89.7 90.6 89.8 88.0 

Southern. 80.5 82.7 82.9 77.6 

X 2 -ratio. 22.6 21.0 19.7 22.3 

Pseudo-collapsed 3-point function 

Northern. >99.9 >99.9 >99.9 99.7 

Southern. 35.1 34.9 35.8 31.4 

X 2 -ratio. 98.8 98.5 98.5 98.4 

Equilateral 3-point function 

Northern. 98.6 98.6 98.8 98.4 

Southern. 45.7 45.7 47.8 42.6 

X 2 -ratio. 86.6 86.7 86.6 86.7 

Rhombic 4-point function 

Northern. 99.7 99.7 99.7 99.6 

Southern. 22.8 22.5 23.2 20.1 

X 2 -ratio. 97.3 97.1 97.2 97.0 


Table 16. Probabilities of obtaining values for the x 2 statistic 
and ratio of x 2 of the N -point functions shown in Fig. 18 for the 
Planck fiducial ACDM model at least as large as the observed 
values of the statistic for the SMICA map on hemispheres de¬ 
fined by the ecliptic (first column), Doppler boost (DB, second 
column), and dipolar modulation (DM, third column) reference 
frames. 


Probability [%\ 


Hemisphere Eel. DB DM 

2-point function 

Negative . 88.0 86.9 61.8 

Positive. 77.6 91.1 59.9 

X 2 -ratio. 22.3 5.1 7.7 

Pseudo-collapsed 3-point function 

Negative . 99.7 64.1 95.9 

Positive. 31.4 79.3 48.3 

X 2 -ratio. 98.4 23.3 78.6 

Equilateral 3-point function 

Negative . 98.4 54.8 >99.9 

Positive. 42.6 95.0 78.4 

X 2 -ratio. 86.7 67.7 88.2 

Rhombic 4-point function 

Negative . 99.6 46.4 97.5 

Positive. 20.1 86.3 23.2 

X 2 -ratio. 97.0 57.9 92.5 


The results for the 7V-point correlation functions deter¬ 
mined in the DB and DM reference frames for the SMICA 
map are shown in Fig. 18 and the probabilities are presented 
in Table 16. Note that the positive hemisphere for the eclip¬ 
tic reference frame corresponds to the southern hemisphere 
in the previous table. Whilst the largest asymmetry is seen 
in ecliptic coordinates, a substantial asymmetry is present 
also for the DM direction. This can be explained by the fact 
that the DM direction is more closely aligned with the south 
ecliptic pole (with a separation of around 47°) than the DB 
direction is. For the DB direction we do not find any sig¬ 
nificant asymmetry. The equivalent results for Commander, 
NILC, and SEVEM are consistent with those shown here. 

In conclusion, the correlation functions for the Planck 
2015 temperature data are consistent with the results pre¬ 
sented in PCIS13. Specifically, we observe that the northern 
hemisphere correlation functions are relatively featureless 
(both the 3- and 4-point functions lie very close to zero), 
whereas the southern hemisphere functions exhibit a level 
of structure consistent with Gaussian simulations. 


5.3. Constraints on quadrupolar modulation 

The most natural extension of the class of statistically 
anisotropic models that we have considered previously in¬ 
volves the quadrupolar modulation of an initially statisti¬ 
cally isotropic CMB sky map. No detection of a correspond¬ 
ing quadrupolar power asymmetry is currently claimed. An 
initial BipoSH analysis of the WMAP 7-year data (Ben¬ 
nett et al. 2011) found evidence of corresponding non-zero 
spectra, and A|° +2 , in ecliptic coordinates. However, 
Hanson et al. (2010) demonstrated that the signal could be 
attributed to an incomplete treatment of beam asymmetries 
in the data, and this was subsequently confirmed in Bennett 
et al. (2013). The corresponding analysis of the Planck 2013 
data indicated consistency with statistical isotropy (Planck 
Collaboration XXIII 2014). 


Here, we proceed further and consider the quadrupolar 
modulation of the primordial power spectrum as suggested 
by Ackerman et al. (2007): 


P(k) = P{k) 


1 + ^2M(^) 

M 


(32) 


Given such a spectrum, the CMB temperature field is ex¬ 
pected to exhibit a correlation between a^ m and 
with A = 0, 2. Therefore, the BipoSH coefficients Affi and 
Aj ^_ 2 are sensitive to g^M • In the limit of weak anisotropy, 
Kim & Komatsu (2013) proposed an optimal estimator for 
Q2M • 


92M 


l £ ( F_1 )««- £ £ 


m' 


dg2M' 
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(33) 


32M=0 


where a is the CMB data vector in harmonic space and C 
is its covariance matrix, and 
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(34) 


52M=0 


Here, {(C~ 1 a*) em (C~ 1 a)e' m f ) g2M=0 is the mean field in 
the absence of the quadrupolar modulation. Observation- 
specific issues such as incomplete sky coverage, inhomoge¬ 
neous noise, and asymmetric beams will result in a non-zero 
mean field, which can be estimated for the Planck data us¬ 
ing simulations. Due to the otherwise prohibitive compu¬ 
tational cost, we adopt a diagonal approximation for the 
inverse of the covariance matrix: 


(C ^)lm,Pm' ~ l/(^€ A Ng) Su'SmmP 


(35) 
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Fig. IT. Difference of the iV-point correlation functions determined from the N S id e = 64 Planck CMB 2015 temperature estimates 
and the corresponding means estimated from 1000 simulations. Results are shown for the 2-point, pseudo-collapsed 3-point (upper 
left and right panels, respectively), equilateral 3-point, and connected rhombic 4-point functions (lower left and right panels, 
respectively). Correlation functions are shown for the analysis performed on northern (blue) and southern (red) hemispheres 
determined in the ecliptic coordinate frame. The solid, dashed, dot-dashed, and dotted lines correspond to the Commander, NILC, 
SEVEM, and SMICA maps, respectively. Note that the lines lie on top of each other. The shaded dark and light grey regions indicate, 
for reference, the 68 % and 95 % confidence regions, respectively, determined from the SMICA simulations. 


where Ct and Ni are the signal and noise power spectra 
respectively. Uncertainties are computed by applying the 
estimator to simulations. 

Table 17 presents results from an analysis of the Planck 
data using the extended common mask, UTA76, and limit¬ 
ing the range of multipoles to 2 < £ < 1200. When including 
data at higher ^-values, simulations show evidence for large 
statistical uncertainties in the recovered g^M values that 
are a consequence of the many holes in the mask related 
to point sources. Therefore, imposing this limit i < 1200 
does not significantly affect the constraining power of the 
analysis. We then estimate the amplitude of the quadrupo- 

lar modulation using the relation g<± = (l/5 \92m\ 2 Y 2 • 

Due to the nature of the estimator, which is necessarily pos¬ 
itive, the estimation is biased. For an unbiased assessment, 
we estimate the mean and standard deviation of g<± from 


simulations. We find no evidence for quadrupolar modula¬ 
tion of the primordial power spectrum. However, the de¬ 
rived limits allow us to impose tight constraints on sta¬ 
tistically anisotropic inflationary models, such as those in¬ 
cluding vector fields during inflation. A companion paper, 
Planck Collaboration XX (2015), contains a more complete 
discussion on the theoretical implications of this constraint. 

5.4. Point-parity asymmetry 

The CMB anisotropy field defined on the sky, T(n), may be 
divided into symmetric, T + (n), and antisymmetric, T _ (n), 
functions with respect to the centre of the sphere, as pre¬ 
viously described in PCIS 13. These functions have even 
and odd parity, and thus correspond to spherical harmon¬ 
ics with even and odd £- modes, respectively. On the very 
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Fig. 18. Difference of the iV-point correlation functions determined from the N S id e = 64 Planck SMICA CMB 2015 temperature 
estimates and the corresponding means estimated from 1000 simulations. Results are shown for the 2-point, pseudo-collapsed 3- 
point (upper left and right panels, respectively), equilateral 3-point, and connected rhombic 4-point functions (lower left and right 
panels, respectively). Correlation functions are shown for the analysis performed on negative (blue) and positive (red) hemispheres 
determined in the ecliptic (solid lines), Doppler boost (DB, dashed lines), and dipole modulation (DM, dot-dashed lines) coordinate 
frames. The shaded dark and light grey regions indicate the 68 % and 95 % confidence regions, respectively. 


Table 17. Constraints on the quadrupolar modulation, determined from the Commander, NILC, SEVEM, and SMICA foreground- 
cleaned maps. The first three columns correspond to the five independent parts of the quadrupolar modulation, which we have 
chosen to present using a complex notation for ^m- The quoted error bars are at the 68% confidence level. The quadrupolar 
modulation amplitude is given in the fourth column, while the mean and standard deviation of $ 2 , estimated from simulations, 
are provided in the fifth column. 


g 2 M X 10 2 g 2 X 10 2 


Method M = 0 M = 1 M = 2 Data Simulation 

Commander. 1.31 ± 1.22 (0.43 ± 0.86) + i (-0.01 ± 0.68) (1.08 ± 0.89) + i (-0.38 ± 0.86) 0.97 1.12 ±0.37 

NILC. 0.88 ±1.21 (0.37 ±0.85) +i ( 0.33 ± 0.67) (0.87 ± 0.88) ± * (-0.26 ± 0.86) 0.77 1.11 ±0.37 

SEVEM. 0.85 ± 1.22 (0.35 ± 0.85) ± i ( 0.34 ± 0.67) (1-00 ± 0.88) ± i (-0.25 ± 0.86) 0.81 1.11 ±0.37 

SMICA. 1.10 ±1.10 (0.46 ±0.81) ±j( 0.26 ± 0.64) (0.93 ± 0.83) ± * (-0-26 ± 0.82) 0.85 1.05 ±0.34 


large scales corresponding to the Sachs-Wolfe plateau of the should be parity neutral with no particular parity prefer- 
temperature power spectrum (2 < £ < 30), the Universe ence exhibited by the CMB fluctuations. However, an odd 
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Fig. 19. The ratio i? TT (Cnax) for Commander (red), NILC (or¬ 
ange), SEVEM (green), and SMICA (blue) determined at N S id e = 
32. The shaded grey regions indicate the distribution of the 
statistic derived from the SMICA MC simulations, with the dark, 
lighter, and light grey bands corresponding to the 1, 2, and 3 a 
confidence levels. 

point-parity preference has previously been observed in the 
WMAP data releases (Land & Magueijo 2005a, b; Kim & 
Naselsky 2010a, b; Gruppuso et al. 2011) and the Planck 
2013 results. Here, we investigate the parity asymmetry in 
the 2015 temperature maps at A S id e = 32. We consider the 
following estimator: 

dTT /a \ _ -P+Tmax) /oa ^ 

-K (/max) — jjTT £ \> O”) 


Fig. 20. Lower-tail probability of the point-parity estimator 
for Commander (red), NILC (orange), SEVEM, (green), and SMICA 
(blue). 


occurs at i = 28 corresponding to a value of 0.2 % for NILC, 
SEVEM, and SMICA, and 0.3% for Commander. 4 

As a first attempt to quantify any a posteriori effects in 
the significance levels, we consider how many MC simula¬ 
tions appear in the lower tail of the MC distribution with 
a probability equal to, or lower than, 0.2 %, for at least one 
4ax value over a specific range. For £ max in the range 3-50, 
the total number of simulated maps with this property is 
less than 20 over 1000 MC maps, implying that, even con¬ 
sidering the LEE, an odd-parity preference is observed with 
a lower-tail probability of less than 2 %. 


where D + (i max ) and D-(£ max ) are given by 

+,■ 


dI t _ 


^ t 

Hot f—2 f 

k , — ^ I'Lrric 


+ 1) ^TT 


(37) 


is the total number of even (+) or odd (—) multipoles 
included in the sum up to f m ax, and Dj T is the temperature 
angular power spectrum computed using a QML estimator 
(Gruppuso et al. 2011). The factor in Eq. (37) 

effectively flattens the spectrum across the Grange of the 
Sachs-Wolfe plateau (up to i = 50) in a ACDM model. 

Figure 19 presents the ratio, i? TT (Gnax), for the 2015 
component-separated maps, together with the distribution 
determined from the SMICA MC simulations which serves 
as a reference for the expected behaviour of the statis¬ 
tic in a parity-neutral Universe. The distributions for the 
other CMB maps are very similar. The four component- 
separation products are in good agreement, indicating an 
odd-parity preference at very large scales for the multipole 
range considered in this test. 

Figure 20 shows the lower-tail probability for the data as 
compared to simulations as a function of £ max . The results 
are in good agreement with those in PCIS13. The cleaned 
CMB maps yield generally consistent profiles which signify 
an anomalous odd-parity preference in the multipole region 
4ax = 20-30. The minimum in the lower-tail probability 


5.5. Mirror-parity asymmetry 

For the Planck 2013 data release, we studied the proper¬ 
ties of the temperature data at a resolution of A S id e = 16 
under reflection with respect to a plane to search for mir¬ 
ror symmetries. Such a symmetry might be connected to 
non-trivial topologies (Starobinsky 1993; Stevens et al. 
1993; de Oliveira-Costa et al. 1996). In Planck Collab¬ 
oration XXIII (2014), we reported evidence for an anti¬ 
symmetry plane, with a perpendicular direction given by 
(/, b ) = (264°, —17°), However, the probability of the results 
was slightly dependent on the method of foreground clean¬ 
ing, with a p -value ranging from 0.5 % for Commander-Ruler 
to 8.9% for SMICA. The same direction was also found in 
the WMAP 7-year data (Finelli et al. 2012), and is close 
to that determined for the dipole modulation in the Planck 
2013 data release (PCIS13), suggesting possible connections 
between the two directional anomalies. 

We now proceed to reanalyse the status of mirror sym¬ 
metries using the Planck 2015 full mission temperature data 
at both A^ide = 16 and A si d e = 32. In order to avoid pos¬ 
sible bias introduced by the use of the Galactic mask 5 the 
results are derived from the full-sky Commander, NILC, and 

4 In the case where we would like to test the probability of 
finding a Universe with either odd or even parity preference, the 
probability would be higher by a factor of about two. 

5 The Galactic mask induces a preferred direction in the analy¬ 
sis of the MC simulation ensemble, which affects the significance 
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SMICA maps described in Sect. 2. For SEVEM, a customized 
map is first produced by inpainting about 3 % of the map 
along the Galactic plane using a diffusive inpainting tech¬ 
nique. This is then smoothed to the appropriate lower res¬ 
olutions for further analysis. Following Finelli et al. (2012), 
we consider the estimators in the pixel domain given by: 


S' ± (n i ) = 


1 ^ 


AW £ 


1 (8T fir 

2 ( ^( n i) ± -jr( n k) 


(38) 


where the sum is over all 7V pix HEALPix pixels, (ST/T)(hj) 
is the CMB temperature anisotropy measured at the pixel 
defined by the unit vector hj , and is the opposite direc¬ 
tion with respect to the plane defined by hi, i.e., 

hk — ftj — 2 (hi • hj)hi . (39) 

Note that we expect to be small if the points on opposite 
sides of the mirror are negatives of each other, and S~ to 
be small when they are the same. 

We compute these quantities for each of the 3072 
(12288) directions defined at resolution Wide = 16 (32), 
and allow the j and k indices to run over all of the pixels 
of the low-resolution full-sky maps. We perform the same 
analysis on 1000 FFP8 simulations and store the minimum 
value of S± for each of these to compute probabilities. The 
results are summarized in Table 18 and Fig. 21. 

We confirm that the full mission Planck temperature 
data at Wide = 16 exhibit the most anomalous mirror anti¬ 
symmetry in the direction (/, b) = (264°,—17°), consistent 
with the result from the 2013 nominal mission data, with 
a probability which ranges from 1.6% for SEVEM to 2.9% 
for Commander. This is within 40° of the preferred direction 
identified by the dipole modulation analysis in Sect. 6.2. 
The corresponding results at Wide = 32 yield approxi¬ 
mately the same direction, (Z, fe) = (264°,—16°), with a 
slightly increased probability, ranging from 0.8 % for SEVEM 
to 1.9% for Commander. 

We also note that the CMB pattern exhibits a mir¬ 
ror symmetry in the direction (Z, 6) = (260°, 48°), consis¬ 
tent with that found in the WMAP 7-year data (Finelli 
et al. 2012), and close to that identified by the solar dipole 
(Planck Collaboration VIII 2015). However, the significance 
of the symmetry pattern is less than in the antisymmetric 
case. 

This extension of the analysis to higher resolution than 
in our previous work shows that the antisymmetry prop¬ 
erty does not seem to be confined to the largest angular 
scales, although we have not attempted to correct for any 
a posteriori choices made in the analysis. The detailed con¬ 
nection of this antisymmetry property to the low-variance 
and hemispherical asymmetry observed on these scales re¬ 
mains an open issue. 


5.6. Local peak statistics 

Local extrema or peaks, as introduced in Sect. 4.5.3, can 
be employed to search for localized anomalies on the CMB 
sky by examining how their statistical properties vary in 
patches as a function of location. 

Initially, we consider a further test for asymmetry by 
examining the differences in the peak distribution when 

of the results determined from the data. See Ben-David V Kovetz 
(2014) for a discussion. 


Table 18. The lower-tail probability for the S± statistics of the 
component-separated maps at Wide = 16 and Wide = 32. 



Probability 

Direction 

Estimator 

[%] 

M) [°] 

Wide = 16 


Commander 


min(5' + ).... 

2.9 

(264.4,-17.0) 

min^ - ) . . . 

12.0 

(260.4, 48.1) 


NILC 


min(S+).... 

2.3 

(264.4,-17.0) 

min^ - ) . . . 

16.8 

SEVEM 

(260.4, 48.1) 

min(S + ). . . . 

1.6 

(264.4,-17.0) 

min^ - ) . . . 

13.5 

SMICA 

(260.4, 48.1) 

min(5' + ).... 

2.7 

(264.4,-17.0) 

min^ - ) . . . 

19.1 

(260.4, 48.1) 

Wide = 32 


Commander 


min(S+). . . . 

1.9 

(264.4,-15.7) 

min(*S _ ) . . . 

10.0 

(265.3, 46.2) 


NILC 


min(S + ). . . . 

1.2 

(264.4,-15.7) 

min(*S _ ) . . . 

10.3 

SEVEM 

(265.3, 46.2) 

min(5 + ).... 

0.8 

(264.4,-15.7) 

min^ - ) . . . 

11.1 

(265.3, 46.2) 


SMICA 


min(5' + ).... 

1.7 

(264.4,-15.7) 

min^ - ) . . . 

11.6 

(265.3, 46.2) 


divided according to orientation with respect to a previ¬ 
ously specified asymmetry direction. In particular, we se¬ 
lect the peaks both in a disc of radius 70° centred on 
(Z, fe) = (225°,—18°) (the positive direction of the dipole 
defined in Sect. 6.2 for SMICA) and in the corresponding 
antipodal disc, then construct the empirical peak height 
CDFs to be compared with the full-sky median FFP8 dis¬ 
tribution, as shown in Fig. 22. For maps filtered with a 
40' FWHM GAUSS filter the distribution of the peaks 
for the positive-direction disc is in general agreement with 
the full sky result, while that for the negative-direct ion is 
marginally different. Moreover, this pattern of behaviour 
is seen over a number of filtering scales, both for the KS 
deviation from the median full-sky simulated CDFs, and 
the spread of extremal values when comparing positive and 
negative regions. We also find that the properties of the 
negative disc affect the p -value results for a full sky KS test 
on data filtered with an SSG84 filter of 500' FWHM, as 
seen in Sect. 4.5.3. 

We can then extend the analysis for the 40' GAUSS- 
filtered data by considering the variation in the peak sta¬ 
tistical properties for a set of discs, each of which is centred 
on a pixel defined at Wide = 256. The simplest statistics 
to consider are the peak number counts. We therefore con¬ 
sider discs of 30° diameter and compute the peak counts 
for each disc. These are then compared to the correspond¬ 
ing peak count CDFs determined from simulations, and the 
upper- and lower-tail probabilities are assigned by counting 
the number of simulations above and below the observed 
counts at the same location. These quantities can then be 
visualized in the form of Wide = 256 sky maps. The derived 
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Fig. 21. Histograms of the *S + (top panel) and S~ (bottom 
panel) statistic. The vertical lines show the minimum value for 
the estimator computed at iV S ide = 32 for Commander (red), NILC 
(orange), SEVEM (green), and SMICA (blue) maps. The grey area 
shows the same quantity computed from 1000 simulated SMICA 
maps. 


—log 10 (UTP) maps for each component-separation method 
are shown in Fig. 23. While we find that the total counts of 
peaks for the sky coverage defined by the common mask is 
consistent with simulations, significant regional variation is 
seen. Indeed, the p-value for certain disc locations drops to 
0.1% (i.e., the sky counts exceed anything seen in simula¬ 
tions). However, one needs to account for the a posteriori se¬ 
lection of significant regions in the determination of the true 
significance. It should also be noted that regional variations 
of the UTP are seen at similar levels when inspecting the 
peak-count statistics maps derived for randomly selected 
realizations of the simulations. Moreover, the significance 
of such peak-counting anomalies is degraded with larger 
disc diameters, and becomes insignificant for counts on the 
full sky. Thus, no significant anomalies can be claimed for 
the peak-count statistics of the Planck data. 

A powerful non-parametric test of statistical isotropy is 
provided by the two-sample KS-deviation between the full 
sky empirical peak height CDF F n (x) (see Eq. 28) and an 
empirical peak height CDF F n > ( x ) derived from a subsam¬ 
ple of the distribution, again defined by the peaks within 
discs of 30° diameter as defined above. The two-sample KS- 
deviation 

/ Tin! 

Knn' = \ -; - - sup \F n t (x) - F n (x)\ (40) 

V n + n' x 
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Fig. 22. KS-deviation of the peak distribution for 70° radius 
discs centred on the positive and negative asymmetry directions 
determined from the SMICA CMB temperature map in Sect. 6.2. 
From top to bottom, the plots correspond to maps filtered with 
a GAUSS kernel of 40' FWHM, an SSG84 filter of 500' FWHM, 
and an SSG84 filter of 800' FWHM, respectively. 


for a partial sky region shares samples between the two 
CDFs, and can be calculated extremely efficiently using 
rank statistics according to 


1 nn f 

max 

i 

r'(i) - 1 

r(i) — 1 

V n + n' 

n' — 1 

n — 1 


where r and r' denote the ranks of a value with index i in 
the full set of n and restricted set of n' samples, respectively. 
Maps of the upper tail probability are then determined by 
comparison with the equivalent quantities computed from 
simulations; —log 10 (UTP) maps are shown in Fig. 24. The 
majority of the selected locations are consistent with the 
full-sky distribution, thus indicating the statistical isotropy 
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of the Planck maps. The most prominent feature in each of 
the local KS-deviation maps appears south of the Galactic 
centre and may be associated with a cold region crossing 
the Galactic plane. However, as with the peak counts, it 
cannot be interpreted as statistically anomalous. 

5.7. The Cold Spot 

Since its discovery in the WMAP first-year data (Vielva 
et al. 2004), the Cold Spot, centred at Galactic coordi¬ 
nates (/, b ) = (210°, —57°) has been one of the most exten¬ 
sively studied large-scale CMB anomalies. In the 2013 re¬ 
lease (Planck Collaboration XXIII 2014), Planck confirmed 
the apparently anomalous nature of this feature in temper¬ 
ature, in terms of the area of the SMHW coefficients on 
angular scales of « 10° on the sky; the 2015 release has 
also confirmed this feature (see Sects. 4.5.2 and 4.5.3). The 
CMB temperature anisotropies around the Cold Spot as ob¬ 
served by Planck are shown in the top panel of Fig. 25. The 
peak merger tree within the Cold Spot region is presented 
in the lower panel of the figure and provides a multiscale 
view of its structure (see Sect. 4.5.4 for details). 

The robustness of the detection of the anomalies dis¬ 
cussed in this paper is a non-trivial issue. For the particu¬ 
lar case of the Cold Spot, this has been reviewed by Vielva 
(2010), and addressed in detail by Cruz et al. (2006), paying 
specific attention to the impact of a posteriori choices. In 
particular, the latter study focused on the original test that 
indicated the presence of this feature on the sky, confirming 
a significance between 1 % and 2 %. An alternative analysis 
of the significance based on two statistical tests with differ¬ 
ent levels of conservativeness was made by McEwen et al. 
(2005), providing values of 0.1% and 4.7%, respectively. 
The statistical significance of the Cold Spot was questioned 
by Zhang & Huterer (2010) who found a low significance 
after performing a study based on different kernels. As dis¬ 
cussed in more detail by Vielva (2010), this result can also 
be interpreted as evidence that not all kernels are neces¬ 
sarily suitable for the detection of arbitrary non-Gaussian 
features. 

The possibility that the Cold Spot arises from instru¬ 
mental systematics (Vielva et al. 2004) or foreground resid¬ 
uals (Liu & Zhang 2005; Cruz et al. 2006) has been largely 
rejected. However, several non-standard physical mecha¬ 
nisms have been proposed as possible explanations. These 
include the gravitational effect produced by a collapsing 
cosmic texture (Cruz et al. 2007), the linear and nonlin¬ 
ear ISW effect caused by a void in the large-scale structure 
(e.g., Tomita 2005; Inoue & Silk 2006; Rudnick et al. 2007; 
Tomita & Inoue 2008; Finelli et al. 2014), a cosmic bub¬ 
ble collision within the eternal inflation framework (Czech 
et al. 2010; Feeney et al. 2011; McEwen et al. 2012), and a 
localized version of the inhomogeneous reheating scenario 
within the inflationary paradigm (Bueno Sanchez 2014). 

Since the other scenarios lack additional evidence, the 
void hypothesis would seem to be the most plausible, de¬ 
pending on the sizes, density contrasts, and profiles as¬ 
sumed in the computations, some of which are not in agree¬ 
ment with either observation (Cruz et al. 2008) or current 
V-body studies (Cai et al. 2010; Watson et al. 2014). How¬ 
ever, Szapudi et al. (2014) have recently detected a large 
void in the WISE-2MASS galaxy catalogue aligned with the 
Cold Spot, with an estimated radius of around 200h -1 Mpc, 
an averaged density contrast of S « — 0.1, and centred on a 


Table 19. Probabilities of obtaining values for the y 2 statistic 
of the angular profiles of the estimators shown in Fig. 26 larger 
than those determined from the data. 

Probability [%] 

Angular profiles Comm. NILC SEVEM SMICA 


Mean. 0.9 0.8 1.0 0.9 

Variance. 40.0 40.0 38.0 42.0 

Skewness. 79.0 82.0 85.0 80.0 

Kurtosis. 75.0 56.0 75.0 77.0 


redshift of z ^ 0.15. Large voids with similar characteristics 
are not unusual in the standard ACDM model (Nadathur 
et al. 2014). In fact, TV-body simulations predict about 20 
such voids in the local Universe (z < 0.5). However, Zibin 
(2014) and Nadathur et al. (2014) indicate that the ex¬ 
pected signal due to the linear and nonlinear ISW effects 
caused by this structure is not large enough to explain the 
temperature decrement associated with the Cold Spot. 

The new Planck data release allows us to further ex¬ 
plore the statistical nature of the Cold Spot. Two previous 
studies (Zhao 2013; Gurzadyan et al. 2014) have claimed 
inconsistencies of the internal properties of the Cold Spot 
with the Gaussian hypothesis, which we re-address here. In 
particular, we consider the small-scale fluctuations within 
a disc-like region of radius ~ 25°. 

Several statistical quantities are computed from the full- 
resolution temperature maps within the Cold Spot region. 
This is divided into a central disc of diameter 1° surrounded 
by a set of 13 concentric annuli with central radii spaced 
in steps of about 2°, thus allowing us to build angular pro¬ 
files for the mean, variance, skewness, and kurtosis. These 
are then compared to specialized CMB realizations, gen¬ 
erated as follows. A set of Gaussian CMB skies is simu¬ 
lated using the FFP8 reference spectrum, and convolved 
with a Gaussian beam of 5' FWHM. As for the FFP8 sim¬ 
ulations themselves, these maps are rescaled, as discussed 
previously. Only those that contain a spot as extreme as 
the Cold Spot at a scale R = 300' in SMHW space are 
retained, and these are rotated such that each simulated 
cold spot is relocated to the actual position of the Cold 
Spot (this ensures that the noise properties are identical 
for both data and simulations). This selection criterion cor¬ 
responds to the characteristic that originally indicated the 
presence of the Cold Spot in the observed sky. As a final 
step, for each remaining CMB simulation a noise realiza¬ 
tion is added, consistent with each component-separation 
method. 

The results are presented in Fig. 26. Focusing on the 
profile of the mean value, it is apparent that the largest de¬ 
viations from the simulations appear on scales around 15°, 
which corresponds to a hot ring structure, as seen in Fig. 25 
and previously discussed in Cayon et al. (2005) and Na¬ 
dathur et al. (2014). Notice that on the smallest scales the 
mean profile is also somewhat deviant with respect to the 
simulations, but this may be connected to selection bias, 
since we are considering CMB simulations containing a spot 
that is at least as cold as the Cold Spot. However, if we con¬ 
sider the distribution of the profiles corresponding to the 
coldest spots instead of the spots as extreme as the Cold 
Spot (removing the bias at the smallest scales) then the 
results do not change substantially (see below). 
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Fig. 23. Map of —log 10 (UTP) for peak counts in the Planck 40' GAUSS-filtered temperature data, where each pixel encodes the 
probability determined for a 30° diameter disc centred on it. 
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Fig. 24. Map of —log 10 (UTP) for the two-sample KS-deviation where each pixel encodes the probability determined for a 30° 
diameter disc centred on it, as computed from the Planck 40' GAUSS-filtered temperature data. 


In order to quantify possible deviations from Gaussian- 
ity, we determine the probability of finding a x 2 value larger 
than that of the data for each statistic, as summarized in 
Table 19. The x 2 value for the data is computed using an 
estimate of the covariance matrix between different radial 
scales determined from the Cold Spot simulations (1000 for 
each component-separation method), and then compared to 
the theoretical x 2 distribution with 13 degrees of freedom. 
The results indicate that the angular profile for the mean 
is poorly described by the simulations, of which less than 
1 % are found to have a higher x 2 than the data (when 
considering the distribution corresponding to the coldest 
spot this probability becomes approximately 2 %). We have 
checked that this deviation is not obviously associated with 
a particular sub-range of angular scales, implying that the 
mean profile is anomalous over the full range considered. 
Conversely, the radial profiles of the higher-order moments 
are compatible with the Gaussian simulations. The latter 
results are then in contradiction with a similar analysis (us¬ 
ing discs instead of rings) by Zhao (2013) for the WMAP 
9-year data. However, it appears that this may be related to 
the criteria applied for the selection of the Gaussian sim¬ 
ulations used to define the null hypothesis. In particular, 
Zhao (2013) used the coldest pixel in real space as a means 
to identify those simulations that should be retained, as op¬ 
posed to the existence of cold spots as extreme as the Cold 
Spot selected in the SMHW coefficient map at R = 300'. 
Since it is not implicit that such a temperature extremum 
is necessarily associated with an extended cold region, par¬ 
ticularly one defined in wavelet space, the simulations used 
by Zhao (2013) did not contain features comparable to the 
nature of the Cold Spot. This explains why the Cold Spot 
seemed to be anomalous when looking at the small-scale 
fluctuations. 

In conclusion, it appears that only the mean tempera¬ 
ture profile of the Cold Spot should be considered anoma¬ 
lous when compared to CMB cold spots that are as statis¬ 


tically extreme. All other measures of its internal structure 
are consistent with expectations. 

As a final remark, we note that the high-pass filtering 
currently applied to the Planck CMB polarization maps 
severely limits the possibility of conducting targeted analy¬ 
ses to discriminate between different possible origins of the 
Cold Spot. For example, no polarization signal would be 
expected in those models producing secondary anisotropies 
due to a gravitational effect, whereas a specific pattern 
might be expected in a bubble collision scenario (Czech 
et al. 2010). Appropriate tests will be pursued in future 
work, once the large-scale polarization data are available. 

6. Dipole modulation and directionality 

In this section, we examine isotropy violation related to 
dipolar asymmetry, various forms of which have been noted 
since the early WMAP releases (Eriksen et al. 2004a). We 
perform a non-exhaustive series of tests in an attempt to 
narrow down the nature of the asymmetry (on the assump¬ 
tion that it is not simply a statistical fluke). First, we will 
briefly describe some similarities and differences between 
the tests that are important for making a proper compari¬ 
son of the results. 

All the tests in this section share in common the fitting 
of a dipole. This is done either by fitting for a dipole explic¬ 
itly in a map of power on the sky (Sects. 6.1 and 6.5), by 
employing Bayesian techniques in pixel space for a specific 
model (Sect. 6.2), or by measuring the coupling of i to £± 1 
modes in the CMB covariance matrix (Sects. 6.3, 6.4, and 
6.6). The differences arise from how the fitted dipoles are 
combined, which determines the specific form of asymmetry 
that the test is sensitive to. 

The tests can be divided into two categories, amplitude- 
based and direction-based. Sections 6.1 to 6.4 are all sen¬ 
sitive to the amplitude of a dipole modulation. Specifically, 
Sect. 6.1 looks for dipole modulation in the pixel-to-pixel 
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Fig. 25. Top: temperature patch centred on the Cold Spot. 
Bottom: Peak merger tree within the Cold Spot region. The 
figure shows a region centred on the Cold Spot location in 
gnomonic projection, with all the peaks in SSG84-filtered maps 
with FWHM ranging from 80 / to 1200 / overlaid on the same 
plot. The size of the coloured circles is proportional to the filter¬ 
ing scale. The colour corresponds to the peak value, normalized 
in units of a for a given filter scale. In both panels the data are 
from the SMICA CMB map at full resolution. 


variance of the data, while Sects. 6.2, 6.3, and 6.4 all search 
for dipole modulation of the angular power spectrum. The 
distinction between these two approaches is mainly one of 
i weighting. 

Sections 6.5 and 6.6 both examine aspects of direction¬ 
ality in the data, where the directions are extracted from 
dipole fits but combined in different ways. Section 6.5 fits 
for dipoles in band power (with similar results for variance) 
and only uses the direction information, while Sect. 6.6 
weights each dipole equally across all scales and uses the 
amplitude information as well. 

The differences between the approaches of these sections 
should be kept in mind when comparing their results. For 
example, although Sects. 6.5 and 6.6 both look for a direc¬ 
tional signal in the data, they are optimized for different 
forms of deviations from statistical isotropy. It is therefore 
unsurprising that they arrive at different results. However, 
the signal found in Sect. 6.5, if not simply a statistical fluke, 
is constrained by the results of Sect. 6.6. 


6.1. Variance asymmetry 

The study of power asymmetry via the local variance of 
the CMB fluctuations was first performed by Akrami et al. 
(2014) for the Planck 2013 and WMAP 9-year tempera¬ 
ture data. The approach was motivated by its conceptual 
and implementational simplicity, its directly intuitive in¬ 
terpretation, and by virtue of being defined in pixel space, 
a useful complementarity to other mostly harmonic-based 
methods. The statistic was computed over patches of dif¬ 
ferent sizes and positions on the sky, and compared with 
the values obtained from statistically isotropic simulations. 
It was found that none of the 1000 available simulations 
had a larger variance asymmetry than that estimated from 
the data. This suggested the presence of asymmetry at a 
statistical significance of at least 3.3 cr, with a preferred di¬ 
rection (Z, b) « (212°, —13°) in good agreement with other 
studies. In this section, we revisit the variance asymmetry 
and report the results of the analysis for the Planck 2015 
temperature maps at full resolution, N si( ^ e = 2048. 

The analysis proceeds as follows. We consider a set of 
discs of various sizes centred on the pixels of a HEALPix 
map defined by a specific N si ^ e value. For each sky map, 
we first remove the monopole and dipole components from 
the masked sky and then compute the variance of the fluc¬ 
tuations on a given disc using only the unmasked pixels. 
This yields a local-variance map at the HEALPix resolution 
of interest. We also estimate the expected average and vari¬ 
ance of the variance on each disc from the simulations and 
then subtract the resulting average variance map from both 
the observed and simulated local-variance maps. Finally, we 
define the amplitude and direction of the asymmetry by fit¬ 
ting a dipole to each of the local-variance maps, where each 
pixel is weighted by the inverse of the variance of the vari¬ 
ances computed from the simulations at that pixel. At all 
stages, we use only the discs for which more than 10 % of the 
area is unmasked, although our results are robust against 
the choice of this value. The computed local-variance am¬ 
plitudes are then used to compare the data with statisti¬ 
cally isotropic simulations. Note that we use only the dipole 
amplitudes of the local-variance maps to measure the sig¬ 
nificance of the asymmetry; the amplitudes of higher mul¬ 
tipoles were shown by Akrami et al. (2014) to be consistent 
with statistically isotropic simulations and we therefore do 
not consider them in the present paper. 

In Akrami et al. (2014), the sensitivity of the method to 
the disc size was assessed using both statistically isotropic 
and anisotropic simulations. The free parameters, i.e., the 
number and size of the discs, were then fixed by these sim¬ 
ulations. It was found that for 3072 patches centred on the 
set of pixels defined at 7V S id e = 16, the simulated asymmetry 
signals were not detected when either very small (raise < 4°) 
or very large (rai SC > 16°) discs were used. 

The former effect is due to a combination of the low 
number of pixels per disc and an insufficient number of 
discs to cover the entire sky when lV si ae = 16 reference 
grids are used. However, it has recently been shown by Ad- 
hikari (2015) that using a larger number of small discs (by 
increasing N S [d e to 32, 64, 128, and 256, depending on the 
disc size) in order to cover the entire sky allows the local- 
variance method to detect the large-scale anomalous asym¬ 
metry as well as the Doppler boost signal from the Planck 
2013 data, at a significance of > 3.3cr. Fantaye (2014) has 
demonstrated that the Doppler boost signal can be detected 
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Fig. 26. From left to right: the mean, variance, skewness, and kurtosis angular profiles computed for rings at radii 0 centred on 
the Cold Spot position for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The expected value obtained from 
the simulations is denoted by the black dashed line and the grey regions represent the 1 a and 2 a intervals. 


at a similar level of significance using needlet bandpass fil¬ 
tering of the data, even with large discs, when simulations 
are deboosted. Here, in contrast to the 2013 analysis, we 
use maps which contain Doppler boosting, for both simula¬ 
tions and data, and therefore we do not detect any Doppler 
boost signal when using a large number of small discs. 

The low observed significance levels when large discs 
are used is due to the cosmic variance associated with the 
largest-scale modes. Motivated by the analysis of Fantaye 
(2014), and in order to address this issue, we also perform 
analyses using a Butterworth high-pass filter, 


m 


(V4) 4 

1 + 4 / 4 ) 4 ’ 


(42) 


centred at multipoles £q = 5, 10, 15, 20, and 30. In addition, 
the filtering of low multipoles allows us to establish the 
contribution of such modes to any detected asymmetry. 

Here, based on the analysis of Akrami et al. (2014), we 
restrict our analysis to those disc sizes for which 3072 discs, 
corresponding to an N si ^ e = 16 map, cover the entire sky, 
i.e., to the range 4°-90°. Consistent results can be obtained 
by choosing other values of N si( ^ e for a given disc size pro¬ 
vided that the entire sky is covered by the discs. Here, for 
simplicity, we work with the same Aside (=16) for all disc 
sizes. 

Our results for the measured amplitude of the variance 
asymmetry, compared to the values from the simulations, 
as well as the corresponding dipole directions, are shown in 
Fig. 27. The p -values are given for different disc sizes and 
in terms of the number of simulations with local-variance 
dipole amplitudes greater than the ones measured from the 
data. Note that since the discs with different sizes used in 
our analysis are correlated, the significance levels are also 
correlated. For this reason we choose to show the p-values 
as a function of disc size instead of combining them into a 
single number. Moreover, it should be noted that the signif¬ 
icance values we present here do not incorporate any correc¬ 
tions to account for the choice of parameters adopted during 
method calibration, specifically the dipole amplitudes and 
directions for the anisotropic simulations that were used to 
fix the range of disc sizes and number of patches. 

It can be seen from the upper panel of Fig. 27 that for 
the unfiltered map the significance of the power asymmetry 
drops quickly when we increase the disc size to radii greater 


than 16°. This is no longer the case, however, when the low¬ 
est multipoles are filtered out. For example, when the filter 
scale is set to = 5, i.e., when the very low multipoles 
which are affected most by cosmic variance are suppressed, 
the variance asymmetry is detected at the 3 a level for all 
disc sizes, as shown in Fig. 27. Table 20 presents the p- 
values of the variance asymmetry using 8° discs and for 
various values of Our results show that variance asym¬ 
metry is detected with a remarkable significance for all disc 
sizes when very low multipoles are filtered out. In addition, 
the variance asymmetry amplitude slowly decreases with 
increasing as seen in the upper panel of Fig. 28. For 
> 20, the dipole amplitude becomes too small and we 
find no significant variance asymmetry. It is interesting to 
note, however, that the dipole directions found for large 
are closely aligned with those found for < 20. 

Table 20. p -values for the variance asymmetry measured by 8° 
discs for the four component-separated temperature maps and 
different high-pass filter scales. The values represent the fraction 
of simulations with local-variance dipole amplitudes larger than 
those inferred from the data. 

p-value [%\ 


4 Comm. NILC SEVEM SMICA 

Unfiltered . . 0.1 0.1 0.1 0.1 

5. <0.1 <0.1 <0.1 <0.1 

10 . < 0.1 < 0.1 < 0.1 < 0.1 

15. 0.1 0.0 0.1 <0.1 

20. 0.4 <0.1 0.3 0.2 

30 . 1.8 0.8 1.8 1.7 


The lower panel of Fig. 27 shows the dipole directions 
we find using different disc sizes and different filter scales for 
SMICA. The dipole directions for the Commander, NILC, and 
SEVEM component-separated maps are very similar to those 
shown. The asymmetry directions found here are consistent 
with those determined by other analyses in this paper. 

In the upper panel of Fig. 28, we show the local- 
variance dipole amplitudes for the 8° discs as a function 
of the central multipole of the high-pass filter, In the 
lower panel of the same figure we show, as an example, 
the mean-subtracted and inverse-variance-weighted local- 
variance map using 8° discs for the Commander component- 
separation method. The pixels of the map are given in terms 
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Fig. 27. Upper panel: p- values for variance asymmetry measured 
as the number of simulations with local-variance dipole ampli¬ 
tudes larger than those inferred from the data, as a function of 
disc radius for the four component-separated maps, Commander 
(red), NILC (orange), SEVEM (green), and SMICA (blue), and for 
unfiltered and high-pass-filtered cases. For the filtered case, 
the Commander curve is covered by the SMICA curve for small 
( 7 *disk < 8) disks, and by the SEVEM curve for large disks 
(^disk > 8). Lower panel: local-variance dipole directions for the 
SMICA map. The colours, as indicated by the colourbar, corre¬ 
spond to different values of the high-pass filter central multipole 
4. The size of a marker disc corresponds, from small to large, 
to the size of the disc used in the analysis, namely 4°, 12°, 20°, 
and 70°. The dipole directions from the Commander, NILC, and 
SEVEM component-separation methods are consistent with the 
case shown here. The \ow-l and WMAP-9 directions are identi¬ 
cal to those in Fig. 35. 


of the lower- and upper-tail probabilities of the values from 
the data compared to the values from the simulations. The 
maps for NILC, SEVEM, and SMICA are very similar. The nu¬ 
merical values of the local-variance dipole amplitudes and 
directions for the Commander method are given in Table 21; 
the values for NILC, SEVEM, and SMICA methods are similar. 

6.2. Dipole modulation: pixel-based likelihood 

In PCIS13 we presented an analysis of the apparent 
anisotropic distribution of large-scale power in the Planck 
2013 temperature data within the parametric framework 
defined by Gordon (2007) and Hoftuft et al. (2009), who 
introduced an explicit dipole modulation field to model po¬ 


rn 

O 
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Fig. 28. Upper panel: local-variance dipole amplitude for 8° 
discs as a function of the central multipole of the high-pass fil¬ 
ter, 4, for the four component-separation methods, Commander 
(red), NILC (orange), SEVEM (green), and SMICA (blue). The grey 
regions, from dark to light, correspond, respectively, to lcr, 2 a, 
and 3 a percentiles from the 1000 FFP8 simulations processed 
by the Commander method. Lower panel: mean-subtracted and 
inverse-variance-weighted local-variance map for the 8° discs 
and for the Commander component-separation method; each pixel 
is given in terms of the lower- and upper-tail probability of the 
measured value on that pixel compared to the values from the 
simulations. The pixels in grey correspond to the centres of the 
8° discs on which the number of unmasked pixels in the full 
resolution map is lower than our threshold. The black curve 
superposed on the map indicates the boundary of the oppos¬ 
ing hemispheres along the asymmetry axis. It is clear that the 
largest fraction of >95 % outliers (red pixels) lie on the pos¬ 
itive amplitude hemisphere of the local variance dipole, while 
the <5% outliers (blue pixels) are on the opposite hemisphere. 
The corresponding maps for NILC, SEVEM, and SMICA are very 
similar to the one shown here. 


tential hemispherical power asymmetry. The following is a 
direct update of that analysis using the Planck 2015 CMB 
data at 7V S ide = 32, retaining the the 2013 common mask 
to explicitly test for consistency with the earlier study. All 
results are found to be in excellent agreement. In the fol¬ 
lowing, we therefore only consider a smoothing scale of 5° 
FWHM as a representative example. This is the highest 
angular resolution accessible for an 7V si de = 32 map. 
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Table 21. Local-variance dipole amplitudes and directions. 
All values quoted here are for 8 ° discs. This table is for the 
Commander component-separation method, but the results are 
similar for the other methods. 


Direction 

4 A a (Z,6)[°] 


Unfiltered . 0.052 ± 0.016 (210,-26) 

5. 0.046 ±0.014 (208,-24) 

10 . 0.040 ±0.014 (199,-16) 

15 . 0.038 ±0.012 (206,-16) 

20 . 0.028 ± 0.010 (202,-18) 

30 . 0.025 ±0.010 (199,-19) 


a A = 2 (Api anc k — ( Affp8 ))5 where Api an ck and Affps are the 
local-variance dipole amplitudes of the data and the FFP 8 sim¬ 
ulations, respectively. The quoted errors are the dispersion of 
the simulation amplitudes. Assuming a pure dipole modulation 
model, A to first order would correspond to the modulation am¬ 
plitude. 

Table 22. Summary of dipole modulation results at a smoothing 
scale of 5° for all Planck 2015 CMB temperature solutions, as 
derived by the brute-force likelihood given by Eq. 43. 


Method 2013 2015 


Dipole modulation amplitude, a 

Commander. 0.078 ± 0.021 0.066 ± 0.021 

NILC. 0.069 ±0.021 0.061 ± 0.022 

SEVEM. 0.066 ± 0.021 0.065 ± 0.021 

SMICA. 0.065 ±0.021 0.066 ± 0.021 

Dipole modulation direction, (/, b) [°] 

Commander. (227, -15) ± 19 (230, -16) ± 24 

NILC. (226,-16) ± 22 (228,-19) ±29 

SEVEM. (227,-16) ±24 (226,-17) ±25 

SMICA. (226,-17) ± 24 (225,-18) ±24 

Power spectrum amplitude, q 

Commander. • • • 0.961 ± 0.025 

NILC. ••• 0.954 ±0.024 

SEVEM. • • • 0.966 ± 0.025 

SMICA. • • • 0.960 ± 0.025 

Power spectrum tilt, n 

Commander. • • • 0.082 ± 0.043 

NILC. ••• 0.077 ±0.043 

SEVEM. ••• 0.077 ±0.043 

SMICA. • • • 0.081 ± 0.043 



Modulation amplitude, a 



Fig. 29. Top : Marginal constraints on the dipole modulation 
amplitude, as derived from Planck 2015 temperature observa¬ 
tions at a smoothing scale of 5° FWHM for Commander (red), 
NILC (orange), SEVEM (green), and SMICA (blue). The plot corre¬ 
sponds directly to Fig. 32 of Planck Collaboration XXIII (2014). 
The Commander, SEVEM, and SMICA posteriors coincide almost 
perfectly both internally, and with the corresponding SMICA 2013 
posterior, shown as a dashed black line. Bottom : Corresponding 
marginal two-dimensional constraints on the \ow-£ power spec¬ 
trum amplitude and tilt, (q, n), defined relative to the best-fit 


Recall first the basic data model adopted in the dipole 
modulation approach: rather than assuming the CMB sky 
to be a statistically isotropic Gaussian field, we allow for 
an additional dipole modulation, resulting in a data model 
of the form d = BMs ± n. where My = (1 ± ap • hi)8ij is 
an offset dipole field multiplying an intrinsically isotropic 
signal s with a dipole of amplitude a pointing towards 
some preferred direction p. B denotes convolution with 
an instrumental beam, and n denotes instrumental noise. 
Additionally, we model the power spectrum of the un¬ 
derlying statistically isotropic field in terms of a two- 
parameter amplitude-tilt model of the form Ct(q,ri) = 
q (^/30) n C ACDM , where (7 ACDM is the best-fit Planck 2015 
ACDM spectrum (Planck Collaboration XI 2015). The two 
parameters q and n can accommodate a deficit in power 
at low £ as compared to the best-fit cosmology that would 
otherwise create a tension with the underlying statistically 


Planck 2015 ACDM model. 


isotropic model and result in the analysis measuring a com¬ 
bination of both asymmetry and power mismatch. 

In the absence of any dipole modulation, a = 0, the 
total data covariance matrix is given by C = BSi SO B T ± N, 
where Si so is the standard statistically isotropic CMB co- 
variance matrix given by the power spectrum, C^, N is the 
noise covariance matrix, and the corresponding likelihood 
is given by the usual expression for a multivariate Gaus¬ 
sian distribution. With dipole modulation, this generalizes 
straightforwardly to C = BMSi SO M T B T ± N, with the likeli¬ 
hood given by 


C(a,p, q , n) oc 


exp [-id*(BMSM T B T + N) -1 d] 
y|BMSM T B T "TN| 


(43) 
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Table 23. Amplitude (A) and direction of the \ow-£ dipole mod¬ 
ulation signal determined from the QML analysis for the range 
£ G [2,64]. The errors are calculated from the cosmic variance 
expected for statistically isotropic CMB realizations. 


Direction 

Method A (l,b) [°] 

Commander. 0.063±g;g?| (213, -26) ± 28 

NILC. 0.0641°;°^ (209, -25) ± 28 

SEVEM. 0.063^o'oi3 (211, -25) ± 28 

SMICA. 0.062 to'™ (213, -26) ± 28 


Figure 29 and Table 22 summarize this five-dimensional 
likelihood in terms of marginal parameters for each of the 
four Planck CMB maps, as evaluated over the common 
mask using the multi-dimensional grid-based Snake algo¬ 
rithm (Mikkelsen et al. 2013). All results correspond to a 
smoothing scale of 5° FWHM, the highest resolution sup¬ 
ported by an N S [^ e = 32 HEALPix grid, but, as in 2013, we 
consider all smoothing scales between 5° and 10° FWHM, 
reaching similar conclusions in each case: The dipole mod¬ 
ulation results derived from the Planck 2015 temperature 
maps are essentially identical to the 2013 results, with im¬ 
proved internal consistency between the four CMB maps 
due to better mitigation of systematic errors. The best- 
fit dipole modulation amplitude at 5° FWHM is 6-7% 
whilst the low-i? power spectrum has an approximately 3- 
5 % lower amplitude compared to the best-fit ACDM pre¬ 
diction. These results are fully consistent with expectations 
given that the Planck 2013 sky maps were already cosmic- 
variance-limited on these angular scales, and the 2015 maps 
differ from the 2013 maps at the level of only a few mi- 
crokelvin (Planck Collaboration IX 2015). 

6.3. Dipole modulation: QML analysis 

In this section we use the quadratic maximum likelihood 
(QML) estimator introduced in Moss et al. (2011) and de¬ 
scribed in Appendix C to assess the level of dipole modu¬ 
lation in our estimates of the CMB sky at A^ si d e = 2048. 
The specific implementation is essentially identical to that 
used in Hanson & Lewis (2009), Planck Collaboration XVII 
(2014), and Planck Collaboration XXVII (2014), and ex¬ 
ploits the fact that dipole modulation of any cosmological 
parameter is equivalent to coupling of £ to £±1 modes in the 
CMB covariance matrix to leading order (see Appendix C). 
Planck Collaboration XX (2015) presents an alternate anal¬ 
ysis for a specific isocurvature model. 

Since we are interested in dipole modulation there are 
three independent estimators. For our particular approach, 
these are a real-valued m = 0 and a complex-valued m = 1 
estimator, and take the form 

_ 6 m — (F/ m Ti + lm)) 

fio (t + mFe+1 

(44) 

X _ 6 ^em^Cee+iBem(Tl m Tg+i m +i — (T)* m Ti+i m +i)) 

1_ /ii J2 e 6Cj e+1 (£ + l)F e F w 

(45) 


ration XVII (2014), where the approximate filter functions 
are also specified. We define 6Cu +1 = dCt/dX+dCt+i/dX, 
where X is the parameter modulated, and Ag m and are 
numerical coefficients (details can be found in Appendix C). 
The factor fi m corrects the normalization for errors intro¬ 
duced by masking: 

f lm mJdSlY? m (n)M(n), (46) 

where M(Q) is the mask. Finally, we correct the direc¬ 
tion for the effects of inhomogeneous noise which is not ac¬ 
counted for in the filtering process, by weighting the X m by 
the inverse of the variance derived from filtered and mean- 
field corrected simulations. 

The physics is readily accessible in this estimator: the 
^-dependence in modulation determined by the parameter 
X is expressed in the 8Cu +1 factor, and the relevant scales 
appear directly in the limits of the sum. We consider the 
estimator over the range £ m [ n = 2 < £ < £ max . The modu¬ 
lation amplitude and direction are then given by 



(47) 




(48) 

(49) 


It is worth re-emphasizing that the quantities A, 6 , and 
are all dependent on the £ range considered. 

As a consequence of the central limit theorem, for suf¬ 
ficiently large £ max the Vs are Gaussian-distributed with 
mean zero, so that the amplitude parameter has a Maxwell- 
Boltzmann distribution. We fit to this distribution for 
4ax > 10 when computing the p- value, so as not to be 
influenced by Poisson noise in the tails of the empirical dis¬ 
tribution (and we have determined that this is a good fit 
to the simulations by applying a KS test). For the case of 
scalar amplitude modulation (i.e., X = A s ), and f? m i n = 2, 
the cosmic-variance-limited expectation for the modulation 
amplitude from statistically isotropic skies is 


A4A / 48 

A s / V T 4)(^ max — 1) 


(50) 


This is the cosmic variance for a scale-invariant dipole mod¬ 
ulation, and gives a more explicit expression than the 
scaling discussed in Hanson & Lewis (2009). 

The top panel of Fig. 30 presents results for the p- value 
of the fitted modulation amplitude as a function of ^ max . 
Note that there are several peaks, at £ « 40 and £ ~ 67 (the 
focus of most attention in the literature), and £ « 240. The 
latter peak, while not previously emphasized, is also present 
in the WMAP results (see Fig. 15 in Bennett et al. 2011). 
It is also interesting to note that a modulation amplitude is 
observed at ^ max ~ 800 that is somewhat lower than what 
one would typically expect for a statistically isotropic sky. 
However, the significance is not at the level of the excess 
dipole modulation at low £ and will not be discussed fur¬ 
ther. The dip at ^ max ~ 67, with a p-value of 0.9-1.0%, 
corresponds to the well-known low-i? dipole modulation. 6 


Here Tg m are C-inverse filtered data and Fg = 6 Actually only SEVEM and SMICA achieve their minimum at 

We adopt the inverse-variance filter from Planck Collabo- £ max = 67, whereas NILC and Commander achieve theirs at 
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Fig. 30. Probability determined from the QML analysis for a 
Monte Carlo simulation to have a larger dipole modulation am¬ 
plitude than the Commander (red), NILC (orange), SEVEM (green), 
or SMICA (blue) data sets, with (top panel) £ m in = 2 or (bottom 
panel) £ m in = 100. No significant modulation is found once the 
\ow-£ signal is removed. We emphasize that the statistic here is 
cumulative and apparent trends in the curves can be misleading. 


could have occurred over any £ range. Since we are looking 
for a large-scale phenomenon, we assume that the analysis 
should include the corresponding low-^ modes and start at 
£ = 2. In order to correct for a posteriori effects we then 
adopt the following scheme. 

1. We calculate the modulation of each simulation on the 
scales 2-£, where £ £ [10,^ max ]. For each simulation we 
find the modulation that gives the smallest probability, 
r] (in the same way that was done for the data). 

2. With the distribution of 77 s given by the simulations we 
then compare this to the data. That is, we calculate 
the probability that one would find oneself in a Hubble 
patch with a modulation amplitude up to £ £ [ 10 , £ max ] 
that is as significant as (or more significant than) the 
modulation in the real data. 



Fig. 31. Probability determined from the QML analysis for ob¬ 
taining a dipole modulation amplitude at least as anomalous as 
the Commander (red), NILC (orange), SEVEM (green), and SMICA 
(blue) data sets, for the range £ £ [10,£ m ax]. The vertical line 
corresponds to -£ max = 132 which was used as the search limit 
in Bennett et al. (2011). The probability grows approximately 
logarithmically with -£ max . This means that the adopted proba¬ 
bility to exceed is fortunately not very sensitive to -£ max , and for 
any reasonable choice is above 10 %. 


Table 23 presents the corresponding dipole modulation pa¬ 
rameters, which are seen to be consistent with previous 
studies. Note that the mean amplitude expected for a set 
of statistically isotropic simulations at this £ max is 2.9% 
(in close agreement with the expected value due to cosmic 
variance, Eq. 50). 

We have therefore determined a phenomenological sig¬ 
nature of modulation for £ = 2-67 with a p-value of 0.9- 
1.0%. If such a signal had been predicted by a specific 
model, then we could claim a significance of about 3 a. How¬ 
ever, in the absence of such an a priori model, we can assess 
how often we might find a 3 a effect by chance, given that it 


4ax = 14 and 240, respectively. Such scatter is expected when 
searching over a large number of possible £ ranges. The recon¬ 
structed amplitudes for each component-separation method are 
well within the error budgets of the estimator. 


If 4ax = 132 (as chosen by Bennett et al. 2011), the 
probability of achieving a modulation as large as the Planck 
data in this range is higher than 10% (see Fig. 31). This is 
in agreement with the findings of the WMAP team (which 
found 10 % and 13 % in the same Grange, using two differ¬ 
ent masks). Here, we do not quote a specific PTE for the 
dipole modulation since it depends on the choice of both 
4 ax (albeit not so sensitively) and Anin (which we have 
decided not to marginalize over). However, it appears to 
be the case that the dipole modulation that we observe is 
quite unremarkable. That is, Gaussian fluctuations in a sta¬ 
tistically isotropic Universe will reasonably often result in 
a dipole modulation with a comparable level of significance 
to that presented here. 

Beyond this, evidence for dipole modulation is found at 
£ « 200-300, with a smaller dip at £ « 500. Given that 
the dipole modulation estimator is a cumulative quantity, 
it is possible that these features are statistically enhanced 


Article number, page 35 of 66 









A&A proofs: manuscript no. planck_2015_iands_final 


by the usual low-^ signal. To test this we analyse the dipole 
modulation as a function of f m ax again, with the restriction 
^min = 100 applied in order to completely remove any low-£ 
influence. The outcome is presented in Fig. 30 (bottom). It 
is clear that even before introducing posterior corrections 
no significant modulation is found, indicating that the p- 
values of the features at i > 100 were indeed exaggerated 
by the \ow-£ modulation. 


6.4. Bipolar spherical harmonics 

In the absence of the assumption of statistical isotropy, the 
CMB two-point correlation function C(hi, 712 ) ^ C(hi -f^) 
can be most generally expanded in the bipolar spherical 
harmonic (BipoSH) basis representation as follows: 

C(ni,n 2 )= ^2 A^iYe^hi) ®Y e2 (n 2 )} L M ■ (51) 

LMt x t 2 

The BipoSH basis functions, {Y^hi) (8) Yi 2 (ti 2 )}lm are 
tensor products of ordinary spherical harmonic functions, 
and the corresponding expansion coefficients are termed 
BipoSH coefficients (Hajian & Souradeep 2003; Hajian & 
Souradeep 2006). The BipoSH basis provides a complete 
representation of any form of statistical isotropy violation 
with the key advantage of separating the angular scale- 
dependence of the signal in spherical harmonic multipoles, 
t, from the nature of the violation indexed in the bipolar 
multipole space by L. Consequently, it is possible to simul¬ 
taneously determine that such a signal is dipolar (L = 1), 
quadrupolar (L = 2), octopolar (L = 3), and so on, in na¬ 
ture and that the power is restricted to specific ranges of 
angular scales. 

The estimation of BipoSH coefficients from CMB maps 
is a natural generalization of the more routinely undertaken 
estimation of the angular power spectrum Ci . To allow a di¬ 
rect connection to the angular power, we further introduce 
a set of BipoSH spectra at every bipolar harmonic moment, 
(L, M), labelled by a difference index d, defined as follows: 


a lm 
A M+d 


aLM 

^U+d 


n L 


(0 < d < L ), 


(52) 


where C'f 1 !^ t m2 are the Clebsch-Gordon coefficients and 

for brevity the notation n^ 2 ...g n = nlLi V + !)• Bi¬ 
poSH spectra, clearly, are then simply a generalized set of 
CMB angular power spectra, with the standard CMB an¬ 
gular power spectrum Ct — A®® being one of them. 7 While 
A®® quantifies the properties of the statistically isotropic 
part of the CMB fluctuations, the additional BipoSH co¬ 
efficients quantify the statistically anisotropic part of the 
CMB two-point correlation function. 

Thus BipoSH provides a mathematically complete de¬ 
scription of all possible violations of statistical isotropy in 

7 The BipoSH spectra, as defined in Eq. (52), restrict us to 
working with only even-parity BipoSH coefficients (L-\-d is even) 
due to the vanishing of Cjj$_|_ d0 otherwise. While most known 
isotropy-violating phenomena like weak lensing, Doppler boost, 
non-circular beams, etc., can only produce even-parity BipoSH 
spectra, measurement of odd-parity BipoSH spectra can be used 
to test for systematic effects, or to search for the signatures of 
exotic effects such as the lensing of CMB photons by tensor 
metric perturbations. 


a Gaussian CMB sky map. It is then always possible to 
translate any specific model for such a signal into the lan¬ 
guage of BipoSH and provide a common approach for the 
multiple specialized tests that have been implemented pre¬ 
viously in this paper and elsewhere. However, improving on 
the analysis of the 2013 Planck data, a new formalism is 
developed in order to reliably analyse a masked sky, as con¬ 
cisely described in Appendix D. Aluri et al. (2015) provides 
a more detailed description of the approach and includes an 
explicit demonstration of its validity using simulations. 

Initially, we revisit the simple phenomenological model 
of dipole modulation of the CMB sky from Sect. 6.2, 

T(h) = To(n) (1 + M(h)) , (53) 

where T(h) represents the modulated CMB sky, Tb(n) is 
the underlying (statistically isotropic) random CMB sky, 
and M(h) is a dipolar field. The BipoSH coefficients re¬ 
sulting from such a modulation are given by 


1 M 

‘tt+1 

= A\^ +1 + mi M G \ e+1 ; 

(54) 

tl 

Ct + Ct+ 1 1 (2£ + 1) (2£ + 3) rl0 

(55) 

'tt+1 

V 3 Oo(^+i)o • 


Here A\#^ 1 corresponds to the BipoSH coefficients of the 
unknown, but statistically isotropic, unmodulated CMB 
field, m im are the spherical harmonic coefficients of the 
modulation field, and Ct is the best-fit CMB angular power 
spectrum. 

The BipoSH representation further enables an estimate 
of the modulation field to be made over specific angular 
scales by windowing regions in multipole space in the sum 
over multipoles £ in Eq. (55). This additional information is 
important for identifying the origin of the isotropy-breaking 
signal, which could be either cosmological or due to system¬ 
atic artefacts. 

We perform the analysis for the N si( i e = 2048 compo¬ 
nent separated CMB maps with an apodized version of the 
common mask at that resolution and reconstruct the mod¬ 
ulation signal in independent bins of width A£ = 64 up to 
4ax = 512. The application of the common mask intro¬ 
duces a mean field bias in the BipoSH coefficients derived 
from the data. This bias is estimated from the FFP8 sim¬ 
ulations and subtracted from the derived coefficients. The 
process of masking induces a coupling between the modu¬ 
lation field and the mask that results in a modification of 
the spectral shape of the modulation signal by the mod¬ 
ified shape function (MSF) (see Appendix D for details). 
Further, the covariance of the bias-subtracted BipoSH co¬ 
efficients is not easy to derive analytically in this case. To 
overcome this problem, we consider the diagonal approxi¬ 
mation to the covariance matrix and estimate it from sim¬ 
ulations. 

The results presented in the top panel of Fig. 32 in¬ 
dicate that the dipole modulation signal is most signifi¬ 
cant in the lowest multipole window £ G [2,64]. Note that 
the power in the dipole modulation field mi = (|ran| 2 + 
|mio| 2 + |mi_i| 2 )/3 is related to the dipole amplitude by 
A = 1.5 y^mi/ir. The best-fit amplitude (A) and direc¬ 
tion corresponding to the reconstructed dipole modulation 
field from this lowest multipole bin is quoted in Table 24 
for each component-separation method. Also shown are 
the corresponding results for the cleaned frequency maps 
SEVEM-100, SEVEM-143, and SEVEM-217. As expected for 
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Fig. 32. Top: Measured dipole modulation (L = 1) power in 
non-overlapping CMB multipole bins for Commander (red), NILC 
(orange), SEVEM (green), and SMICA (blue) as determined from 
a BipoSH analysis of the data. The power in the dipole of the 
modulation field is a y 2 -distributed variable with 3 degrees of 
freedom. The shaded regions in the plot depict, in dark-grey, 
grey, and light-grey respectively, the 1, 2, and 3cr equivalent 
intervals of the distribution function derived from simulations, 
while the solid black line denotes its median. Significant power in 
the dipole modulation is seen to be limited to £ — 2-64 and does 
not extend to higher multipoles. Bottom: Dipole modulation di¬ 
rection as determined from the SMICA map. The directions found 
from the other component separation maps are consistent with 
this analysis. The coloured circles denote the central value of 
the multipole bin used in the analysis, as specified in the colour 
bar. The low-^ and WMAP-9 directions are identical to those in 
Fig. 35. 


signals with a cosmological origin, no evidence for frequency 
dependence is seen. 

Since the amplitude of the dipole modulation field is 
consistent with zero within 2 a for all of the higher Gbins 
considered, it is plausible that the simple modulation model 
in Eq. (53) is inadequate to describe the features seen in the 
BipoSH spectra and should minimally allow for the ampli¬ 
tude, A(£), of the dipole to depend on CMB multipole, £. 
Although this may appear to be a more complex model, it 
does not necessarily lack motivation. It is readily conceiv¬ 
able that physical mechanisms that cause a dipolar modu¬ 
lation of the random CMB sky would be scale-dependent 
and possibly significant only at low wavenumbers. It is also 


Table 24. Amplitude ( A ) and direction of the dipole modula¬ 
tion in Galactic coordinates as estimated for the multipole range 
£ E [2,64] using a BipoSH analysis. The measured values of the 
dipole amplitude and direction are consistent for all maps. 


Direction 

Method A (l,b) [°] 

Commander.. 0.067 d= 0.023 (230,-18)4=31 

NILC. 0.069 ±0.022 (228,-17) ±30 

SEVEM. 0.067 ±0.023 (230,-17) ±31 

SMICA. 0.069 ±0.022 (228,-18) ±30 

SEVEM-100 . . 0.070 ± 0.023 (231, -19) ± 30 

SEVEM-143 . . 0.068 ± 0.023 (230, -17) ± 31 

SEVEM-217 . . 0.069 ± 0.023 (229, -20) ± 31 


intriguing to note that, although in most cases the ampli¬ 
tude of the modulation dipole is seen at low significance, the 
directions in the first four bins, £32 E [2, 64], £qq E [65,128], 
^160 £ [129,192], and ^224 £ [193,256], are seen to be clus¬ 
tered together, as shown in the bottom panel of Fig. 32. 
Note that the lower significance of the modulation for the 
multipole bins at £ > 64 results in larger errors for their re¬ 
spective directions than the value quoted for the £ E [2, 64] 
bin recorded in Table 24. 

We extend our analysis to carry out the dipole modu¬ 
lation reconstruction in cumulative bins up to £ max = 512, 
making cumulative increments in the multipole in steps of 
A£ = 64. The results of this analysis are summarized in 
Fig. 33. 

As noted previously, as a consequence of our motion 
with respect to the CMB rest frame, the observed CMB 
map is expected to be statistically anisotropic, as has been 
demonstrated in Planck Collaboration XXVII (2014) and 
Appendix B. Reassuringly, in PCIS13 it was established 
that such a signal would not contaminate a dipole modula¬ 
tion signal up to £ max ~ 700. We now confirm the Doppler 
boost signal using the BipoSH methodology. 

An equivalent description of the Doppler boost in terms 


of BipoSH 

coefficients is given by 



aim 

= A% a +/3i M G} lia , 


(56) 

GU 

= {MgU]*MgU 

f}x 



/(2£ 1 + l)(2£ 2 + l) 
V 127 T 

r>10 

Wi(M 2 0 ^ 

(57) 


= [C 4 + Cp 2 \ , 


(58) 

[GUd 

= [cy + Ci 2 ] 


(59) 


+ — Cg 2 ] [r 1 (r 1 -+ 

- 1 ) — ^ 2(^2 

+ l)]/ 2 , 


where /3im = f dnYi M (n) ft • n, f3 = v/c denotes the pe¬ 
culiar velocity of our local rest frame with respect to the 
CMB, and b v is the frequency-dependent boost factor, as 
discussed in more detail in Planck Collaboration XXVII 
(2014). 

Since the Doppler boost signal has a frequency de¬ 
pendence, we perform our analysis on the SEVEM-100, 
SEVEM-143, and SEVEM-217 maps at V side = 2048, and 
adopt values of b v = 1.51,1.96, and 3.07, respectively. A 
minimum variance estimator for f3i m, as discussed in Ap¬ 
pendix D, is adopted with the shape function Gf l£ replaced 


Article number, page 37 of 66 




























A&A proofs: manuscript no. planck_2015_iands_final 



100 200 300 400 500 

£ 



64 128 192 256 320 384 448 512 

P 

^max 


Fig. 33. Top: Measured dipole modulation power in cumula¬ 
tive CMB multipole bins for Commander (red), NILC (orange), 
SEVEM (green), and SMICA (blue) as determined from a BipoSH 
analysis of the data.. Colour coding as in Fig. 32. Note that 
the measurements in cumulative bins indicate a power in excess 
of 2 a up to multipole -£ max ~ 320. The value on the horizon¬ 
tal axis denotes the maximum multipole used in the analysis, 
with Cnin = 2. Bottom: Modulation dipole direction as recov¬ 
ered from the SMICA map. The directions found from the other 
component-separation maps are consistent with these directions. 
The colour-coded points represent the directions recovered for 
the specific ^ max used in the analysis, with £ min = 2. The \ow-£ 
and WMAP-9 directions are identical to those in Fig. 35. 

by the corresponding Doppler boost term given in Eq. (56). 
Corresponding unboosted CMB simulations were also used, 
in particular to correct for the mean field bias. However, we 
use a set of Doppler-boosted simulations in order to esti¬ 
mate the error on the reconstructed Doppler boost vector. 

Since it is expected that the low multipole modes of the 
i spectrum are contaminated by the dipolar signal re¬ 
ported previously, in order to monitor the impact of this 
anomalous signal on the Doppler reconstruction we imple¬ 
ment a cumulative analysis using multipoles with a vary¬ 
ing ^ min from 2 to 640 in increments of A^ min = 128 and 
a fixed f m ax = 1024. 8 The recovered Doppler amplitudes 
from the three SEVEM frequency cleaned maps as a func¬ 


8 We fix f max = 1024 since at higher £ values the mismatch 
between the data and simulation power spectra becomes more 


tion of f? m in are shown in the top panel of Fig. 34, while 
the lower panel indicates the corresponding direction ft in 
Galactic coordinates determined from the SEVEM-217 data. 
Table 25 records the best-fit amplitudes and directions for 
£ e [640,1024]. 




2 128 256 384 512 640 


min 

Fig. 34. Top: Amplitude \/3\ of the Doppler boost from the 
SEVEM-100, SEVEM-143, and SEVEM-217 maps for different mul¬ 
tipole bins determined using a BipoSH analysis. The maximum 
multipole of each bin is fixed at -£ max = 1024, while £ m in is in¬ 
cremented from £ = 2 to £ = 640 in steps of A£ = 128. The 
dashed line corresponds to the actual dipole boost amplitude, 
|/3| = 1.23 x 10 3 . Bottom: Doppler boost direction j3 measured 
in Galactic coordinates from SEVEM-217. The coloured circles 
denote £ m in used in the analysis, while -£ max = 1024 is held 
fixed. The \ow-£ and WMAP-9 directions are identical to those 
in Fig. 35. 


6.5. Angular clustering of the power distribution 

In the Planck 2013 data release we reported a possible devi¬ 
ation from statistical isotropy in the multipole range £ = 2- 
600, thus confirming earlier findings based on the WMAP 
data (Hansen et al. 2009; Axelsson et al. 2013). This claim 
of asymmetry extending to higher multipoles was made only 
on the basis of the alignment of preferred directions as de¬ 
termined from maps of the power distribution on the sky 


important and is a concern for the bias subtraction applied when 
reconstructing the Doppler boost signal. 


Article number, page 38 of 66 

































































Planck Collaboration: Isotropy and statistics of the CMB 


Table 25. The Doppler boost amplitude (|/3|) and direction 
in Galactic coordinates derived over the multipole range i E 
[640,1024] as evaluated from a BipoSH analysis. The errors are 
estimated from an identical analysis of a set of 1000 Doppler 
boosted simulations for each frequency. 


Method 

Direction 

101 xio- 3 M)[°] 

SEVEM-100 . . . 

SEVEM-143 . . . 

SEVEM-217 . . . 

. 1.24 ±0.66 (277,40) ±50 
. 1.35 ±0.56 (264,39) ±39 
. 1.28 ±0.45 (257,42) ±32 


for specific multipole ranges. In particular, it was found 
that the directions of the dipoles fitted to such maps in the 
multipole range i = 2-600 were significantly more aligned 
than in simulations. In addition, we showed that the ratio of 
the power spectra in the two opposite hemispheres defined 
by the asymmetry axis for i = 2-600 was not statistically 
anomalous (as later confirmed over the extended multipole 
range i = 2-2000 by Quartin & Notari 2014). 

Here, we test for the alignment in the Planck 2015 data 
set. We adopt the approach for the estimation of the dipole 
alignment that was described in detail in PCIS13, a brief 
summary of which follows. 

1. Local power spectra are estimated from the data at 
lV S ide = 2048 for 12 patches of the sky corresponding 
to the TVside = 1 HEALPix base pixels. Only those high- 
resolution pixels surviving the application of the com¬ 
mon mask are included in the analysis. 9 As a conse¬ 
quence of this masking, when patches based on HEALPix 
pixels with N S [^ e > 1 are used, the available sky fraction 
for those patches close to the Galactic plane is too small 
for power-spectrum estimation. For most of the analysis, 
we use the cross-spectra determined from half-mission 
data sets. 10 Due to a mismatch between the noise level 
in the data and the simulated maps, the results based 
on auto-spectra are less reliable and also more prone 
to other systematic effects than the cross-spectra. We 
therefore do not consider such results here. The spectra 
are binned over various bin sizes between A £ = 8 and 
M = 32. 

2. For each power spectrum multipole bin, an A si d e = 1 
HEALPix map with the local power distribution is con¬ 
structed. 

3. The best-fit dipole amplitude and direction are esti¬ 
mated from this map using inverse-variance weighting, 
where the variance is determined from the local spec¬ 
tra computed from the simulations. We do not compute 
error bars for the direction, but expect this to be ac¬ 
counted for in part by the use of equivalently treated 
simulations in the clustering analysis. 

4. A measure of the alignment of the different multipole 
blocks is then constructed. In PCIS13, we considered 

9 Departing from the analysis in PCIS13, we do not use an 
apodized version of the common mask. Simulations indicate that 
the error on the power spectrum for those multipoles in the range 
300 to 500 where the significance is highest is up to 20 % larger 
in this case, with the corresponding error on preferred direction 
being typically 8 % larger. 

10 Note that simulated half-mission noise maps were generated 
by adjusting the properties of the existing 1000 (10 000 in the 
case of SMICA) noise simulations appropriately, thus explaining 
why only 500 (5000) simulations are used in this analysis. 


the mean angle between all possible pairs of dipole di¬ 
rections up to a given ^ max . Here, for greater consistency 
with Sect 6.6, we use the mean of the cosine of the 
angles, rather than of the angles themselves, between 
all pairs of dipoles. This effectively corresponds to the 
Rayleigh statistic (RS) introduced formally in Sect. 6.6, 
and we will refer to it as such, although it differs by ig¬ 
noring all amplitude information. Clearly, smaller values 
of the RS correspond to less clustering. 

5. The clustering as a function of ^ max is then assessed us¬ 
ing p -values determined as follows. We first construct 
the RS using all multipoles up to f m ax- The p -value is 
then given by the fraction of simulations with a higher 
RS than for the data for this ^ max - A small p -value there¬ 
fore means that there are few simulations that exhibit 
as strong clustering as the data. Note that the p-values 
are highly correlated as the RS is a cumulative function 
of Q 

Ui l ma X - 

6. We then define two measures of significance. To achieve 
this, it is necessary to reduce the 1499 different p -values 
determined for ^ max E [2,1500] to a single measure of 
clustering. We do this in two different ways, using the 
mean of these p- values, and by finding the minimum of 
the p- values, for both the data and for each available 
simulation. We then determine the percentage of sim¬ 
ulations with (i) a lower mean p-value and (ii) a lower 
minimum p-value than the data. Note that these two 
measures of significance take into account different as¬ 
pects of the data. Note further that since the RS is cu¬ 
mulative and the p -values therefore correlated, different 
scales are weighted unequally and a detection in the 
mean and/or minimum p -value may be difficult to in¬ 
terpret and to correct for the multiplicity of tests effect 
(LEE). 

Note that the statistics defined in step 6 above corre¬ 
spond to two choices of what were referred to as “global 
statistics” in PCIS13 in order to assess the degree to which 
the significance of the results depends on a specific choice 
for ^ max . The mean p -value over all available £ max mea¬ 
sures the degree to which clustering is present over large 
multipole ranges independently of whether the clustering is 
strongly focused in one given direction. Clearly the p -values 
for different ^ max are strongly correlated, but if the cluster¬ 
ing is present only over a small multipole range, the RS will 
drop and the corresponding p-values will eventually rise. 
By comparing this value to simulations, we test not only 
whether the dipole alignment in the data is stronger than in 
statistically isotropic random simulations, but also whether 
it is present over larger ranges of multipoles than expected. 
The minimum p- value will give strong detections if there is a 
strong asymmetry over a limited multipole range or weaker 
clustering over larger multipole ranges when the clustering 
is strongly focused in a given direction. 

For Commander, NILC, and SEVEM, only 500 simulations 
are available. However, 5000 simulations are available for 
SMICA, which allows a better estimate of significance to be 
determined when the probabilities obtained are very low. In 
this case, we use half of the 5000 simulations to calibrate the 
statistic (obtain p -values following step 5 above) and the re¬ 
maining half to determine significance levels (compute the 
mean and minimum over these p-values as a function of 
f max following step 6). When using 500 simulations, it is 
necessary to use the same set of simulations to calibrate as 
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well as to obtain probabilities. A related issue with these 
results is that this set of simulations (corresponding to the 
first 500 out of the 5000 available for SMICA) are observed 
to yield higher p-values for the clustering angle due to a 
statistical fluctuation. Another 9 sets of 500 simulations 
that can be obtained from partitioning the 5000 available 
SMICA simulations all result in lower p- values. As a conse¬ 
quence, we observe that results based on the larger number 
of simulations often give lower p- values than when only 500 
simulations are used. 

In Fig. 35 we show the dipole directions of the 15 lowest 
100-multipole bins for the SMICA map. Here, the binning has 
been chosen for visualization purposes; in further analysis 
of the Planck data we use finer ^-intervals. The preferred 
low-^ modulation direction determined in Sect. 6.2 is also 
indicated, along with the WMAP-9 result determined over 
the range £ = 2 to 600 (Axelsson et al. 2013). The observed 
clustering of the dipole directions is similar to that shown 
in figure 27 of PCIS13. Note that differences in masking, 
foreground subtraction, and residual systematic effects will 
displace the direction of a given dipole with respect to the 
previous analysis. Similar behaviour is seen for all of the 
Planck component-separated maps. 

In PCIS13, we calculated the mean angle between all 
possible pairs of dipole directions determined from maps 
of the local power in multipole bins of size A£ = 16. Here 
we test the possible bias arising from such a choice by con¬ 
sidering bin sizes between A£ = 8 and A£ = 32 in steps 
of 2. The lower limit avoids significant bin-to-bin coupling 
in the power spectra for smaller binnings, whilst the upper 
limit excludes cases where there are an insufficient number 
of derived dipoles from which the mean angle can be calcu¬ 
lated, this leading to poor statistics. In addition to showing 
results for each bin size, we also calculate the variance- 
weighted mean of the power spectra over all bin sizes (the 
Ci for a given bin size is weighted by 1/ y/N^ where 7V b 
is the bin size). In this way, we marginalize over bin sizes 
to obtain local power spectra and thereby the RS for each 
single multipole. 

Figure 36 shows the p- values for the different 
component-separated maps, derived as described in step 
5 above. We see that the results based on 500 simulations 
for NILC, SEVEM, and SMICA are in good agreement. The 
Commander results are less consistent, but this may be re¬ 
lated to the fact that component separation was performed 
independently for the half-mission solutions, in contrast to 
the other methods, where component-separation solutions 
were obtained from the full mission data only. For SMICA, we 
also show p-values based on 2500 simulations. These more 
accurate results show lower p- values, and may indicate that 
those determined from only 500 simulations are not suf¬ 
ficiently stable. Note also that for £ < 100 the p-values 
are not consistent with the detection of a \ow-£ asymme¬ 
try/modulation, as seen by other methods in this paper. 
However, for £ < 100, there are very few bins and the vari¬ 
ance of the RS might therefore be too high for this effect 
to be visible. 

In agreement with the conclusions in PCIS 13, a large 
degree of alignment is seen at least to £ max ~ 600. How¬ 
ever, in contrast to the earlier results where the p-values 
started increasing systematically for £ max > 1000, the 
current p- values remain low for £ max > 750. The full 
component-separated maps which have higher resolution 
and sensitivity are used for the current analysis, instead of 


the single-frequency foreground-cleaned map (SEVEM-143) 
used in PCIS13. We note that the results for the updated 
SEVEM-143 map are consistent with the earlier analysis, 
both with and without correction for the Doppler mod¬ 
ulation. Note also that the SMICA results with improved 
statistics (based on 2500 simulations) generally show lower 
p-values than the corresponding results based on 500 sim¬ 
ulations. 

Table 26 presents the fraction of simulations with a 
lower mean/minimum p- value than in the data for a num¬ 
ber of different cases. The table shows probabilities for 
SMICA with different bin sizes (showing only every second 
bin size since these are correlated), as well as for the re¬ 
sults marginalized over bin sizes. We also show results for 
the different component-separated maps, results based on 
half-ring cross-spectra instead of half-mission cross-spectra, 
and results using a different ^-weighting scheme, specifically 
(2^ + l)CV instead of £(£+l )C^, the former being a measure 
of the variance of the temperature fluctuations. The ta¬ 
ble indicates probabilities of approximately 0-2 % for most 
of these cases, although results for the smallest bin size 
show much less significant results. This could be due to 
the strong anticorrelations between adjacent bins found for 
this bin size in those Galactic A^e = 1 patches with very 
small available sky fraction. For the other bin sizes, these 
correlations are much weaker. Note that many of the sig¬ 
nificances based on minimum p- value are only upper limits. 
This is due to the fact that the limited number of simula¬ 
tions in some cases results in the lowest minimum p- value 
being zero. When the minimum p-value in the data is zero, 
we show the percentage of simulations which also have zero 
as the minimum p-value. Clearly this fraction is only an 
upper limit on the real significance. 

In order to further investigate the ^-dependence of the 
asymmetry, we follow two approaches from PCIS 13. Firstly, 
we restrict the analysis to multipoles above a minimum 
multipole £ m { n . Table 26 indicates that clustering at the 
< 1 % significance level is still found when considering only 
those multipoles with ^ min greater than 100. However, when 
this limit is increased to 200, no significant clustering is 
found. We then calculate the RS between pairs of dipoles 
where one dipole is determined from an Grange above a 
certain limiting multipole Gim, and the other dipole below 
this limit. Figure 37 shows the RS as a function of ^ max 
for some selected values of £u m . The £u m = 300 curve (pur¬ 
ple) indicates that dipole directions for £ > 1000 are sig¬ 
nificantly aligned with dipoles for £ < 300. Similarly, the 
£iim = 700 curve (cyan) indicates that the dipole directions 
for £ = 700-1000 are strongly correlated with the dipole 
directions for £ < 700. 

Combining these results, we note that when using only 
multipoles with (i) £ > 200, or (ii) £ < 200, no significant 
clustering is found. The strong clustering significance shown 
to persist to high multipoles in Fig. 36 must therefore be the 
result of clustering of the dipole directions between low and 
high multipoles as supported by Fig. 37. The low p- values 
can be explained by the alignment of dipole directions for 
multipoles extending all the way to £ = 1500 correlated with 
directions for £ < 200. The observed asymmetry is therefore 
not consistent with a model based on dipole modulation or 
power asymmetry located in one specific multipole range 
or for one given direction, but rather as a correlation of 
the dipole directions between £ < 200 and £ > 200. This 
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Fig. 35. Dipole directions for independent 100-multipole bins of the local power spectrum distribution from £ — 2 to 1500 in the 
SMICA map with the common mask applied. We also show the preferred dipolar modulation axis (labelled as “low-A’) derived in 
Sect. 6.2, as well as the total direction for -£ max = 600 determined from WMAP-9 (Axelsson et al. 2013). The average directions 
determined from the two multipole ranges £ E [2,300] and £ E [750,1500] are shown as blue and red rings, respectively. The error 
on the derived direction that results from masking the data is about 60°, with only small variations related to bin size. 
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Fig. 36. Derived p -values for the angular clustering of the power 
distribution as a function of An a x, determined for Commander 
(red), NILC (orange), SEVEM (green), and SMICA (blue), based 
on 500 simulations. For SMICA, the p-values based on 2500 
simulations are also shown (black). The p -values are based on 
the fraction of simulations with a higher RS, determined over 
the Arange up to the given An ax , compared to the data. The 
results shown here have been marginalized over bin sizes in the 
range A£ = 8 to A£ — 32. 


correlation with lower multipoles is found to persist all the 
way to £ max = 1500. 



Fig. 37. Derived p- values for the angular clustering analysis as 
a function of An ax , determined from SMICA based on 2500 sim¬ 
ulations. The p- values are based on the fraction of simulations 
with a higher Rayleigh statistic up to the given £ max than in 
the data. The RS here is calculated over all pairs of dipole di¬ 
rections where one dipole in each pair is computed in the range 
[Aim, An ax ], and the other is determined in the range [2, Aim]. The 
plot shows p- values for Aim = 300 (purple), Aim = 400 (yellow), 
Aim — ^00 (pink), and Aim = 700 (cyan). The results have been 
marginalized over bin sizes in the range A£ — 8 to A£ — 32. 


An advantage of the directional analysis performed here 
is that it focuses on a central issue for tests of deviation from 
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isotropy — whether there is a preferred direction. Indeed, 
Bunn & Scott (2000) noted that the CMB may exhibit a 
pattern that cannot be identified from the power spectrum, 
but which would indicate some non-trivial large-scale struc¬ 
ture. Evidence for the close correlation and alignment of 
directions on different angular scales may present a signa¬ 
ture of broken statistical isotropy, since in the standard 
model, these directions should all be independent random 
variables. In this context, we do not quote a specific direc¬ 
tion for such asymmetry here since our results indicate a 
clustering of angles between different multipoles, but not 
necessarily that all multipoles are clustered about one spe¬ 
cific direction. However, crucially we have shown that the 
measured clustering is driven by the correlations of direc¬ 
tions between higher and lower multipoles. 

Some of the analyses in other sections of the paper fo¬ 
cus on dipolar modulation, a specific model for a dipolar 
power enhancement of the statistically isotropic CMB field 
towards a preferred direction of the sky, and use methods 
optimized for the detection of such a signal. While the re¬ 
sults of Sect. 6.6 show no detection of the clustering of 
directions, there is no clear contradiction with the results 
presented here, since they are based on tests for correla¬ 
tions between different multipoles as expected in the dipolar 
modulation model. The clustering analysis presented here 
is a model-independent test for deviations from statistical 
isotropy which could induce very different correlation struc¬ 
ture. It is therefore sensitive to other forms of asymmetry, 
such as the addition of power in one part of the sky or more 
general phase correlations. 

6.6. Rayleigh statistic: QML analysis 

Results from Sect. 6.5 and in PCIS 13 suggest that, beyond 
a dipole modulation of power on large angular scales, some 
form of directional asymmetry continues to small scales. 
There are also indications from Sect. 6.5 that the direc¬ 
tions of dipolar asymmetry are correlated between large and 
small angular scales. Since the nature of the asymmetry is 
unknown we use the RS, a generic test for directionality 
that makes minimal assumptions about the nature of the 
asymmetry. This statistic has been used both in previous 
CMB studies (Stannard & Coles 2005) and other areas of 
cosmology (Scott 1991). In our context, for a statistically 
isotropic sky this statistic is identical to a three-dimensional 
random walk. The implementation here incorporates all in¬ 
formation pertaining to modulation, not just the direction. 
The approach in this section differs from that of Sect. 6.5 in 
the method of reconstructing power, the choice of binning, 
and the choice of how to weight directions in each bin. An¬ 
other important difference is that Sect. 6.5 only considers 
the direction of dipolar asymmetry and does not take into 
account its amplitude. 

The statistic is cumulative and thus narrowing down the 
specific scales from which a signal may be originating is a 
non-trivial task. However, it is the case that all statistics 
that measure this form of asymmetry (dipole modulation or 
large-scale clustering of power) are in some way cumulative 
and so we will not worry about this issue any further. An¬ 
other disadvantage of this approach is that it will generally 
be less powerful than a test that uses a specific model for 
the directionality. Again, this is a distinction shared when 
one compares any non-parametric versus parametric statis¬ 
tic. 


Table 26. Significance of the angular clustering of the power 
distribution. We indicate the actual mean/min p -value of the 
data, determined from Fig. 36 and written as a fraction of the 
number of simulations used to assess the values, together with 
the percentage of simulations with a lower mean/minimum p- 
value than the data. Unless otherwise specified, the numbers 
are determined from half-mission cross spectra Ct£{£ + 1), for 


all multipoles 
mask. 

in the 

range £ — 

2-1500, 

and for the 

common 


Bin 

Mean 

% 

Min. 

% 

Method 

size 

p -value 

(mean) 

p -value 

(min) 

SMICA . 

8 

261/2500 

1.60 

35/2500 

16.2 

SMICA . 

10 

51/2500 

0.08 

3/2500 

2.36 

SMICA . 

12 

75/2500 

0.20 

1/2500 

0.96 

SMICA . 

14 

83/2500 

0.16 

2/2500 

1.52 

SMICA . 

16 

78/2500 

0.24 

4/2500 

2.00 

SMICA . 

18 

51/2500 

0.04 

1/2500 

0.68 

SMICA . 

20 

21/2500 

<0.04 

1/2500 

0.76 

SMICA . 

22 

60/2500 

0.08 

2/2500 

1.24 

SMICA . 

24 

34/2500 

0.08 

2/2500 

1.00 

SMICA . 

26 

38/2500 

0.08 

1/2500 

0.96 

SMICA . 

28 

42/2500 

0.20 

0/2500 

<0.52 

SMICA . 

30 

27/2500 

0.20 

0/2500 

<0.60 

SMICA . 

32 

21/2500 

0.04 

0/2500 

<0.52 

SMICA . 

marg. 

43/2500 

<0.04 

0/2500 

<1.00 

SMICA a . 

marg. 

48/2500 

<0.04 

1/2500 

1.70 

SMICA b . 

marg. 

47/2500 

<0.04 

0/2500 

<1.16 

SMICA c . 

marg. 

50/2500 

<0.04 

0/2500 

<0.76 

SMICA d . 

marg. 254/2500 

1.52 

34/2500 

20.1 

Comm. 

marg. 

9/500 

<0.20 

0/500 

<2.60 

NILC. 

marg. 

10/500 

<0.20 

0/500 

<3.60 

SEVEM . 

marg. 

13/500 

<0.20 

0/500 

<4.00 

SMICA . 

marg. 

11/500 

<0.20 

0/500 

<3.60 

Comm. b . 

marg. 

11/500 

<0.20 

0/500 

<3.00 

NILC b . 

marg. 

10/500 

<0.20 

0/500 

<3.80 

SEVEM b . 

marg. 

12/500 

<0.20 

0/500 

<3.40 

SMICA b . 

marg. 

11/500 

<0.20 

0/500 

<3.80 

Comm. c . 

marg. 

8/500 

0.20 

0/500 

<4.00 

NILC c . 

marg. 

14/500 

0.20 

1/500 

7.20 

SEVEM C . 

marg. 

17/500 

0.20 

1/500 

8.40 

SMICA c . 

marg. 

15/500 

0.20 

1/500 

7.60 


a Half-ring maps instead of half-mission maps. 
b Ci( 2£ + 1) instead of Cg£(£ + 1). 
c Restricted to multipoles £ > 100. 
d Restricted to multipoles £ > 200. 

The construction of the statistic is as follows. 

1. Beginning with the estimator from Eqs. (44) and (45), 
compute the following binned quantities for the data 
and simulation: 

SC u+1 F e F e+1 (i+l) 1 (60) 

_ 6 B£m{Tg m T£+i rn+1 ~ m+l)) 

fn 4 1) 

(61) 

For each £ this computes the coupling of £ to £ + 1. We 
emphasize that this is a very natural choice of binning 
the estimator, since any parameter that is dipole mod¬ 
ulated will lead to coupling of £ to £ ± 1 modes, albeit 
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with different ^-weightings (below we describe why this 
is not an important issue). 

2. Construct a three-dimensional vector out of the three 
estimators for both the data and the simulations, 11 as 
defined by Eqs. (47-49). 

3. Compute the mean amplitude from simulations and di¬ 
vide all vectors (data and simulations) by this ampli¬ 
tude. This choice ensures that each vector is treated 
equally, since we have no a priori reason to weight some 
scales more than others. 

4. Add this new vector to the previous vector. If this is the 
first time going through this process the previous vector 
is the zero vector. 

5. Repeat with £ £ + 1. Note that the statistics of 

this process are identical to a three dimensional random 
walk. 

Given that a dipole modulation amplitude of roughly 3 a 
significance is known to exist at low £ (before a posteriori 
correction), one would expect a similar level of detection 
of asymmetry to be determined by the RS. Indeed, we find 
that asymmetry is present out to £ « 240. Figure 38 (top) 
presents the p -values derived when the RS is computed as a 
function of ^ max from £ = 2. The minimum p -value obtained 
by the data is 0.1-0.2%, to be compared to the value of 
0.9-1.0% obtained for the dipole modulation amplitude at 
4ax = 67. The direction preferred by the data for ^ max ~ 
240 is (Z, 6) = (208°,—29°), which is approximately 20° 
away from the dipole modulation direction determined to 
£^64. 

We correct for a posteriori statistics using the same pro¬ 
cedure as in Sect. 6.3. Specifically, we count how often sim¬ 
ulations find asymmetry in the range 10 < £ < f m ax that is 
more significant than that found for the data. From Fig. 39 
it is clear that generic asymmetry at the significance level 
found in our CMB sky occurs about 6 % or 8 % of the time 
(depending on the range of £ one decides to search over). 

While the PTE here is not very low, it is nevertheless 
somewhat lower than for the usual dipole modulation test. 
Hence, it seems worth exploring whether any of this sig¬ 
nal comes from higher multipoles. Therefore we compute 
the RS starting at £ m i n = 100, to avoid the influence of 
asymmetry at lower £. The lower panel of Fig. 38 presents 
the corresponding p -values as a function of ^ max . There is 
a striking similarity with the lower panel of Fig. 30. It is 
clear that, even in the absence of a posteriori correction, 
we find no significant asymmetry at larger £. Hence most of 
the signal we are seeing in Fig. 38 (top) is due to the usual 
low-£ asymmetry. 

We would like to stress that the results here are very 
similar to the results of the previous section. For each of 
the statistics used we are simply asking whether there is 
significant coupling of £ with £±1 modes. The details of how 
to optimally combine these couplings for a given £ range 
depends on whether we are talking about dipole modulation 
or directionality (or some other related test, e.g., variance 
asymmetry). These details will change the range of scales 
over which the strongest signal in the data is found. 

11 Note that here we have not specified what 5Cu +1 is (it is 
fully specified by choosing a parameter X to modulate). This 
is because we have decided to weight each £ equally and thus 
any strictly positive choice for 8Cu +i will be equivalent, since 
in step 3 we force the mean length of the vectors at each £ to be 
equal. 




Fig. 38. Rayleigh statistic p -values determined from the QML 
analysis as a function of £ ma x for the Commander (red), NILC 
(orange), SEVEM (green), and SMICA (blue) data sets, with (top 
panel) £ m in = 2 and (bottom panel) £ m in = 100. The general 
pattern of peaks is very similar to that in Fig. 30. We emphasize 
that the statistic here is cumulative and as such trends in the 
curves can be misleading. 

7. Sensitivity of anomalies to enhanced sky 
coverage 

One of the critical aspects in searching for anomalous fea¬ 
tures in sky maps is to ensure that the region being in¬ 
vestigated constitutes a fair and unbiased sample. Since 
many of the claimed anomalies are on large angular scales, 
this implies that minimal masking should be applied to 
the data. However, residual foregrounds then become a 
significant consideration. The masks applied to the four 
component-separated maps studied in the bulk of this pa¬ 
per have been defined at high resolution, and then con¬ 
servatively degraded for lower resolution studies. Such a 
procedure inevitably reduces the sky coverage available for 
analysis, and can be particularly problematic if significant 
structures are aligned by chance with the masked regions. 
Indeed, the WMAP team (Bennett et al. 2011) have drawn 
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Fig. 39. Probability to exceed (PTE) the p -value of the signal 
from the Commander (red), NILC (orange), SEVEM (green), and 
SMICA (blue) data at £ = 230-240 (which is the multipole range 
with the most significant deviation) when searching over a range 
of multipoles up to £ ma x, for the RS determined from the QML 
analysis. Much like the equivalent curve for dipole modulation, 
the PTE appears to grow approximately logarithmically with 

^max- 

attention to several such features in their ILC reconstruc¬ 
tion of the CMB sky, and these are clearly also present in 
the Planck Commander, NILC, SEVEM, and SMICA sky maps. 
A large cold spot is seen near to the Galactic centre, a 
significant fraction of which lies within the common mask 
at any resolution. However, despite its location and visual 
impression, the feature is neither likely to be attributable 
to residual foreground emission, nor is it inconsistent with 
the ACDM model (Gott et al. 2007). In addition, four elon¬ 
gated cold fingers stretching from near the Galactic equator 
to the south Galactic pole are seen, although no equiva¬ 
lent features are evident in the northern sky. Bennett et al. 
(2011) have noted that the alignment of the £ = 2 and £ = 
3 multipoles (Tegmark et al. 2003) seems to be intimately 
connected with these large-scale cool fingers and the inter¬ 
vening warm regions. One of the latter also corresponds to 
the well-known “Bianchi Vllh” main lobe originally found 
in Jaffe et al. (2005). 

Although we would ideally pursue full sky analyses, we 
prefer to remain mindful of the influence of residual fore¬ 
grounds, but still seek to minimize the extent of any mask 
applied for analysis. In this context, and specifically for 
large-angular-scale studies, we consider the properties of an 
additional estimate of the CMB sky, also generated using 
the Commander component separation methodology. In par¬ 
ticular, we note that the Planck low-^ likelihood analysis 
(Planck Collaboration XI 2015) uses the temperature solu¬ 
tion from this study, degraded to a resolution of N s id e = 16. 
The Lkl-Commander map, as we now refer to it, is ini¬ 
tially derived from input data sets (32 bands) at 1° FWHM 
resolution and A^ide = 256. This includes Planck indi¬ 
vidual detector and detector set maps from 30-857 GHz, 
the 9-year WMAP observations between 23 and 94 GHz, 
and the 408 MHz sky survey (Haslam et al. 1982), whereas 
the Commander map described in Planck Collaboration IX 
(2015) includes Planck data alone. It is believed that the 


32-band solution is better (on large angular scales) than the 
Planck-only map, because the larger number of input fre¬ 
quencies allows more detailed foreground modelling, and in 
particular the separation of the low-frequency foregrounds 
into synchrotron, free-free, and spinning dust components. 
An associated confidence mask (hereafter LklT 2 5693) is 
then defined based on a goodness-of-fit measure per pixel, 
corresponding to a rejection of 7.3 % of the pixels on the 
sky. A detailed discussion of these results can be found in 
Planck Collaboration X (2015). 

We now consider the implications of using the 
Lkl-Commander map for studies of several large-angular- 
scale anomalies observed in previous sections, in particu¬ 
lar since the larger sky coverage permitted by this data 
set should constitute a better sample of the Universe. Note 
that, at the resolutions of interest for the following analyses, 
the noise level is negligible (even accounting for the WMAP 
contribution) and should not have significant impact on the 
results. The exact details of the noise contribution to sim¬ 
ulations is therefore unimportant. 

7.1. Variance, skewness, and kurtosis 

We begin by estimating the variance, skewness, and kur¬ 
tosis of the CMB. We apply the unit variance estimator 
to the Lkl-Commander map, and specifically to the version 
used in the low-^ likelihood analysis, which is smoothed to 
440' FWHM at a resolution of N si a e = 16. A corresponding 
\ow-£ mask is generated by a simple degrading of the mask 
at A/gide = 256, then setting those A^ide = 16 pixels with a 
value less than 0.5 to zero and all others to unity. The resul¬ 
tant \ow-£ likelihood mask rejects only 6.4 % of the sky. We 
compare the results for both this mask (also to be referred 
to as LklTi 6 94), and the standard common mask at this res¬ 
olution (UTie58). The results are summarized in Table 27 
and show that, when using the \ow-£ likelihood mask, the 
lower tail probability for the variance is 7.0%. This value 
is higher than the corresponding values for the component 
separated maps as shown in Table 12. In addition the skew¬ 
ness and kurtosis are less consistent with Gaussianity than 
the component separated maps. However, when using the 
standard common mask at N S id e = 16, the lower tail proba¬ 
bility of the variance, skewness, and kurtosis become more 
compatible with those derived earlier. 

There are two possible explanations for this behaviour. 
Either the variance of the CMB in the region close to the 
Galactic plane is intrinsically high, perhaps due to the pres¬ 
ence of the various features noted above, or the presence 
of residual foregrounds increases the variance of the map. 
In order to attempt to distinguish between these options, 
we again apply the unit variance estimator to the stan¬ 
dard component-separated maps 12 , but this time utilising 
the \ow-£ mask. Although the component-separated maps 
are likely to contain some foreground contamination in the 
regions omitted by application of the UTi658 mask, it is 
appropriate to recall that this was constructed in a con- 

12 Note that the SEVEM maps used in this section have been in- 
painted within 3% of the sky towards the Galactic centre us¬ 
ing a simple diffusive inpainting algorithm. This prevents resid¬ 
ual foreground contamination from propagating to neighbour¬ 
ing regions when downgrading the map. The other component- 
separated maps are not pre-processed in this way since some 
form of inpainting of the most contaminated regions was already 
implemented as part of the component separation algorithms. 
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Table 27. Lower-tail probability for the variance, skewness, and 
kurtosis of the Lkl-Commander map. 


Probability [%] 


Mask Variance Skewness Kurtosis 

LklTi 6 94 . 7.0 1.5 94.0 

UTi 6 58 . 0.7 19.9 82.5 


Table 28. Lower-tail probability for the variance, skewness, and 
kurtosis of the Lkl-Commander map compared to the compo¬ 
nent separated maps, obtained using the low-£ likelihood mask 
LklTi 6 94. 


Probability [%] 


Map Variance Skewness Kurtosis 


Lkl-Commander. 7.0 1.5 94.0 

Commander . 7.7 1.9 96.0 

NILC. 9.6 5.0 94.4 

SEVEM. 7.4 4.8 94.3 

SMICA. 7.7 3.7 93.7 

SEVEM-100 . 8.4 0.4 97.9 

SEVEM-143 7.7 3.7 95.5 

SEVEM-217 . 8.3 0.7 95.2 


servative way, and may also mask parts of the sky where 
the level of residual foregrounds can be considered negligi¬ 
ble. In addition, we investigate the cleaned frequency maps 
produced by the SEVEM algorithm in order to test for the 
presence of frequency-dependent residual foregrounds. The 
results of the unit variance estimator analysis are summa¬ 
rized in Table 28. 

All of the component separated maps show an increase 
in the lower tail probability from about 0.5 % when the 
UTi 6 58 mask is applied to roughly 7% for the LklTi 6 94 
mask. The small variations in results for the different maps 
may be attributable to the presence of residual foregrounds 
close to the Galactic plane. However, the increased proba¬ 
bilities can also be explained by the presence of CMB struc¬ 
tures with higher variance within that region which is not 
rejected by the less conservative mask. Indeed, since the 
component-separated maps are affected by different resid¬ 
ual foregrounds, if the source of the changes in probabilities 
is due only to the residual foregrounds, then we would ex¬ 
pect a larger dispersion than what is observed. We also note 
that when we apply the \ow-£ likelihood mask the skewness 
and kurtosis values are shifted towards more extreme val¬ 
ues. This implies that the sky signal is less Gaussian for the 
larger sky fraction, despite the results remaining compati¬ 
ble with the ACDM model assumed for the null tests. Both 
Commander maps are noteworthy in this regard. 

An important issue is whether the changes in the statis¬ 
tics can simply be attributed to differences in the masks. We 
determine how many simulations show an increase in vari¬ 
ance at least as large as that seen for the Lkl-Commander 
map when comparing the values derived for the UTi658 and 
low-^ likelihood masks. Similarly, we determine how many 
simulations have increased skewness or kurtosis values with 
shifts at least as large as observed. When the three statistics 
are considered separately, the fraction of simulations that 
indicate such changes are 7.6%, 4.3%, and 13.9% for the 


variance, skewness, and kurtosis, respectively. Of course, 
such subsets of the simulations also include cases where 
a large shift in the statistic is observed, but the statistic 
would not be considered anomalous for either mask. If we 
also impose the requirement that the simulations have these 
shifts for all three quantities simultaneously, then only 2 
maps from 1000 are found. Of course, such a requirement 
is rather strong, and at this stage we are likely to be ap¬ 
proaching the limits of what can be said based on model- 
independent null tests. Indeed, in order to assess whether 
these results are sensitive to a posteriori choices, we repeat 
the analysis but successively take each simulation as the ref¬ 
erence. Thus, for each simulation the shift in the variance, 
skewness, and kurtosis is computed and then we determine 
how many times we find a case in which two or less of the 
remaining simulations simultaneously show larger shifts for 
the three moments. We find that 48 maps from 1000 satisfy 
these conditions. Given this, it is difficult to draw strong 
conclusions about the significance or otherwise of the mask- 
related changes in variance. 

7.2. N-point correlation functions 

The connection between sky coverage and the observed 
structure of the 2-point correlation function for large an¬ 
gular separations has previously been discussed in the lit¬ 
erature, in particular in connection with the Si/ 2 statistic 
discussed in Sect. 5.2. Bennett et al. (2011) consider that 
the use of a Galactic mask when computing these quanti¬ 
ties is sub-optimal, and note that a full-sky computation 
of the 2-point correlation function from the 7-year WMAP 
ILC map lies within the 95 % confidence region determined 
by simulations of their best-fit ACDM model over all an¬ 
gular separations. However, Copi et al. (2009) suggest that 
the origin of the inconsistencies between the full-sky and 
cut-sky large-scale angular correlations remains unknown, 
and that the observed discrepancies may indicate that the 
Universe is not statistically isotropic on these scales. We 
therefore consider the 7V-point correlation functions, and 
related statistics, of the Lkl-Commander map to contribute 
to this debate. 

We compare results computed for both the 
Lkl-Commander and Commander maps at N s ^ e = 64 

after smoothing to a FWHM of 160'. A mask is con¬ 
structed for the Lkl-Commander map by degrading the 
LklT25693 mask to 7V S id e = 64 and setting all resulting 
pixels with a value less than 0.5 to zero, with the remainder 
set to unity. The LklT 64 92 mask retains 92 % of the sky, to 
be compared to the 67% usable sky coverage allowed by 
the UT6467 common mask at this resolution. 

The results are presented in Fig. 40 where we compare 
the 7V-point functions for the data and the mean values esti¬ 
mated from 1000 Commander simulations. The probabilities 
for obtaining values of the % 2 statistic for the Planck fidu¬ 
cial ACDM model at least as large as the observed values 
are provided in Table 29. For the estimation of the probabil¬ 
ities, we use the same set of 1000 Commander simulations for 
both versions of the Commander data. As noted previously, 
the details of the simulations for such highly smoothed data 
is essentially unimportant. We also provide an analysis of 
the Lkl-Commander map using the common mask to enable 
a direct comparison with the analysis of the Commander 
map. In this latter case, the results for both maps are in 
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Table 29. Probabilities to exceed the observed values of the x 2 
statistics for the Lkl-Commander and Commander maps at iV S ide = 
64. 


Function 

Probability [%] 


Lkl-Comm. a 

Lkl-Comm. b 

Comm. b 

2 -pt. 

84.3 

97.1 

97.2 

Pseudo-coil. 3-pt. 

76.8 

92.1 

92.1 

Equil. 3-pt. 

96.5 

74.0 

74.0 

Rhombic 4-pt. 

94.5 

65.0 

64.6 


a low-f mask, LklT6492 
b common mask, UT6467 


Table 30. Probabilities for obtaining values of the S i/ 2 and Xo 
statistics for the simulations at least as large as the observed 
values of the statistic estimated from the Lkl-Commander and 
Commander maps using the LklT6492 and UT6467 masks, respec¬ 
tively. We also show the corresponding estimation of the global 
p -value for the S(x) statistic. 



Probability [%] 

Statistic 

Lkl-Comm. Comm. 

Sl/2 . 

_ 97.1 99.5 

S(x) (global) . 

_ 90.9 97.7 

xl{0> 60°) . 

_ 95.7 98.1 


excellent agreement. However, the Lkl-Commander map is 
more consistent with simulations when the LklT6492 mask 
is adopted for the 2-point and pseudo-collapsed 3-point 
functions, but less consistent for the equilateral 3-point and 
rhombic 4-point function results. Nevertheless, the results 
are generally in agreement with expectations for a Gaus¬ 
sian, statistically isotropic model of the CMB fluctuations. 

The increased consistency of the 2-point function with 
simulations when analysing a larger sky fraction is consis¬ 
tent with the observations in Copi et al. (2009). We there¬ 
fore quantify this further by determining the statistical 
quantities introduced in Sect. 5.2 for the Lkl-Commander 
map. In particular, we reassess the lack of correlation de¬ 
termined previously for large angular scales. It is evident 
from Table 30 that the results for the Si/ 2 and Xo statistics 
are less anomalous when the low-£ mask is applied. More¬ 
over, the global p- value for the S(x) statistic is substantially 
smaller. 

We also repeat the conventional x 2 analysis but con¬ 
straining the computations to the two separate ranges de¬ 
fined by 6 < 60° and 6 > 60°. The results of these studies 
are shown in Table 31. The analysis for seperation angles 
0 > 60° indicates that the unusually good fit of the observed 
2-point function to the mean 2-point function determined 
for the ACDM model is independent of the mask used in 
the analysis. Conversely, the results for the angles 0 < 60° 
indicate a strong dependence on the mask. It appears that 
the decreased significance of the x 2 statistic for the 2-point 
function of the Lkl-Commander map reported in Table 29 
is related mainly to correlations in the data for separation 
angles smaller than 60°. 

Our results do appear to indicate that computations 
made on larger sky fractions increase the consistency of the 
2-point function with simulations. We therefore also test 
how the hemispherical asymmetry observed previously is 


Table 31. Probabilities for obtaining values of the y 2 statistic 
for the simulations at least as large as the observed values of 
the statistic estimated from the Lkl-Commander and Commander 
maps using the LklT6492 and UT6467 masks, respectively. 



Probability [%] 

Statistic 

Lkl-Comm. Comm. 

x\e < 60°). 

_ 52.9 91.5 

X 2 (0 > 60°). 

_ 96.5 96.8 


Table 32. Probabilities for obtaining values of the x 2 statistic 
and ratio of x 2 of the iV-point functions for the Planck fiducial 
ACDM model at least as large as the observed values of the 
statistic on the northern and southern ecliptic hemispheres esti¬ 
mated from the Lkl-Commander and Commander maps using the 
LklT 6 492 and UT 6 467 masks, respectively. 


Probability [%] 


Hemisphere Lkl-Comm. Comm. 


Northern 

2 -point function 
. 98.7 

89.7 

Southern 

. 29.4 

80.5 

X 2 -ratio . 

. 93.3 

22.6 

Pseudo-collapsed 3-point function 


Northern 

. 99.4 

>99.9 

Southern 

. 14.9 

35.1 

X 2 -ratio . 

. 98.6 

98.8 

Northern 

Equilateral 3-point function 
. 99.9 

98.6 

Southern 

. 47.6 

45.7 

X 2 -ratio . 

. 94.4 

86.6 

Northern 

Rhombic 4-point function 
. 99.8 

99.7 

Southern 

. 24.5 

22.8 

X 2 -ratio . 

. 97.8 

97.3 


affected. The results for the ecliptic frame are presented 
in Table 32. We find that the asymmetry is larger for the 
Lkl-Commander map than for the Commander map in the 
case of the 2-point function, but does not change substan¬ 
tially for the 3-point and 4-point functions. 

7.3. Dipole modulation and directionality 
7.3.1. Variance asymmetry 

Here we apply the local-variance analysis of Sect. 6.1 to the 
Lkl-Commander map and compare the results with those of 
the Commander map. Contrary to the analysis of Sect. 6.1, 
where full-resolution (N si ^ e = 2048) maps were used, here 
the Commander map is downgraded to N s ^ e = 256 in or¬ 
der to consistently compare the results for both maps. The 
simulations used for estimating the significance levels are 
also downgraded to the same resolution, and convolved 
with the corresponding beam function. Otherwise, the pro¬ 
cedure is identical to the one described in Sect. 6.1, e.g., 
the same number of discs has been used to construct the 
local-variance maps. Here we only present the results when 
no high-pass filtering has been applied to the maps; this is 
to avoid confusion as our objective in this section is only to 
compare the general properties of the Lkl-Commander map 
to those of the standard component-separated maps. 

Table 33 summarizes the significance levels measured 
by our variance asymmetry analysis using discs of different 
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Fig. 40. iV-point correlation functions determined from the N S id e = 64 Planck CMB 2015 temperature estimates. Results are shown 
for the 2-point, pseudo-collapsed 3-point (upper left and right panels, respectively), equilateral 3-point, and connected rhombic 
4-point functions (lower left and right panels, respectively). The brown three dot-dashed, purple dashed, and red dot-dashed lines 
correspond to the Lkl-Commander map analysed using the \ow-£ and common masks and the Commander map analysed using the 
common mask, respectively. Note that the dashed and dot-dashed lines lie on top of each other. The black solid line indicates the 
mean for 1000 MC simulations. The shaded dark and light grey regions indicate the 68 % and 95 % confidence regions, respectively, 
estimated using 1000 Commander simulations. See Sect. 4.3 for the definition of the separation angle 6. 


radii, for the Planck 2015 Commander and Lkl-Commander 
temperature maps. The p-values represent the fraction of 
simulations with local-variance dipole amplitudes larger 
than those inferred from the data. We in addition present 
in Table 34 the preferred variance asymmetry directions for 
both maps using 8° discs. 

Our results show consistency between the two maps. 
The small change in the preferred direction is expected from 
the change in the mask, and agrees specifically with the 
directions found by the analysis of the QML dipole mod¬ 
ulation analysis in Sect. 7.3.3. One interesting observation 
is that the large variance asymmetry significance is now 
extended to cases where larger discs are used. Note that 
no high-pass filtering has been applied in the present anal¬ 
ysis, and therefore p-values inferred from the Commander 
map increase with the disc size. As explained in Sect. 6.1, 


the low observed significance levels for larger discs is due 
to the cosmic variance associated with the largest-scale 
modes. The observed increase in the significance levels for 
the Lkl-Commander map is therefore interestingly consis¬ 
tent with this picture; the mask in this case is smaller and 
therefore a larger fraction of the sky is available. This in 
turn provides more data on the largest scales, and there¬ 
fore lowers the impact of the cosmic variance. 

7.3.2. Dipole modulation: pixel-based likelihood 

Table 35 presents constraints on the dipole modulation 
model as derived from the Lkl-Commander map and the 
LklT3293 mask that includes 93% of the sky, updating the 
results from Sect. 6.2 for the Commander map. We find that 
all previously reported results are robust with respect to 
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Table 33. p -values for the variance asymmetry measured by dif¬ 
ferent discs from the Planck 2015 Lkl-Commander and Commander 
temperature solutions using the LklT 25693 and UT25673 masks, 
respectively. The values represent the fraction of simulations 
with local-variance dipole amplitudes larger than those inferred 
from the data. No high-pass filtering has been applied to the 
maps. 


Disc radius [°] 

p-value [%] 

Lkl-Comm. Comm. 

4 . 

0.2 

0.7 

6. 

0.4 

0.3 

8. 

0.3 

0.1 

10. 

0.3 

0.7 

12. 

0.3 

0.8 

14. 

0.5 

1.5 

16. 

0.6 

1.9 

18. 

0.8 

2.7 

20. 

1.1 

3.7 

Table 34. Local-variance dipole directions measured by 8° discs 

for the Planck 2015 Lkl-Commander and Commander temperature 

solutions. 



Method 

M) 

i n 

Lkl-Comm. 

. (225, 

-28) 

Commander . 

. (214, 

-24) 


Table 35. Summary of dipole modulation results at a smoothing 
scale of 5° for the Planck 2015 Lkl-Commander and Commander 
temperature solutions, as derived by the brute-force likelihood 
given by Eq. 43. The former results were derived using the 
LklT 32 93 mask, whereas the latter are those determined pre¬ 
viously in Sect. 6.2. 


Method 

2013 2015 

Dipole modulation amplitude, a 

Lkl-Comm. 

0.059 ±0.020 

Commander . . . . 

. 0.078 ±0.021 0.066 ±0.021 

Dipole modulation direction, (Z, b) [°] 

Lkl-Comm. 

(223, -17) ± 23 

Commander . . . . 

. (227,-15) ± 19 (230,-16) ±24 

Power spectrum amplitude, q 

Lkl-Comm. 

0.970 ±0.025 

Commander . . . . 

0.961 ±0.025 

Power spectrum tilt, n 

Lkl-Comm. 

0.068 ± 0.045 

Commander . . . . 

0.082 ± 0.043 


data selection and sky coverage. In particular, the best-fit 
dipole modulation amplitude at 5° FWHM is 5.9% in the 
Lkl-Commander map, and is thus stable to within about 
0.3cr when increasing the sky fraction from 78% to 93%. 
Likewise, the marginal low-£ power spectrum amplitude, 
q , shifts upward by 0.4 cr, and the power spectrum tilt, n, 
downward by 0.3 cr, for the same sky fraction increases. 

To assess the statistical significance of these shifts, 
we compare with Gaussian statistics, creating two Gaus¬ 
sian random vectors with 78 and 93 elements, respectively, 
where the first 78 elements of the latter vector are iden¬ 
tical to the first vector. From these, we compute the dif¬ 
ference between the two means, after normalizing each so 
that their individual errors in the mean are unity. Repeat- 


Table 36. Summary of the dipole modulation results for 
the range £ G [2,64] determined from the Planck 2015 
Lkl-Commander and Commander temperature solutions, as de¬ 
rived by the QML estimator defined in Sect. 6.3 using the 
LklT25693 and UT78 masks, respectively. 


Direction 

Method A (Z, b ) [°] 

Lkl-Comm. 0.058(227, -28) ± 26 

Commander. 0.063±g;2?f (213, -26) ± 28 


ing this simple calculation 10 5 times, we find that 48 % of 
all Gaussian realizations observe shifts larger than 0.3 cr, 
and 34% observe shifts larger than 0.4 cr. Thus, the param¬ 
eter differences due to the different data selection and sky 
fractions reported above are consistent with expectations 
from random Gaussian statistics. 

7.3.3. Dipole modulation: QML analysis 

We also repeat the QML dipole modulation analysis of 
Sect. 6.3 for the Lkl-Commander map and corresponding 
mask. Table 36 summarizes the results of the low-^ dipole 
modulation for the Lkl-Commander temperature solution, 
compared with the Commander map. 

The best-fit modulation amplitude for Lkl-Commander 
is 5.8 % and the small 0.5 % shift from the Commander best- 
fit amplitude corresponds to a decrease of approximately 
0.4cr. These results mirror very closely the results found 
above for the pixel-based likelihood approach to dipole 
modulation, as expected, and the observed shifts are per¬ 
fectly consistent with those expected from the change in 
the mask. 

7.3.4. Bipolar spherical harmonics 

We next perform a dipole modulation analysis on the 
Lkl-Commander temperature map using the BipoSH for¬ 
malism from Sect. 6.4. The dipole modulation amplitude 
inferred from the analysis is smaller that that deduced from 
analysing the Commander map as seen in Table 37. However, 
it should be noted that the probability for simulations to 
yield a dipole modulation amplitude equal to or greater 
than the amplitude inferred from data is 0.4%, which is 
smaller by a factor of approximately 2.4 as compared to the 
p -value inferred from analysis on Commander. The reduction 
in the dipole amplitude and the enhanced significance can 
both be attributed to the reduced power bias which is a 
result of the increased sky coverage. 

Table 37. Amplitude ( A ) and direction of the dipole mod¬ 
ulation in Galactic coordinates as estimated for the multipole 
range £ E [2,64] using the BipoSH analysis on Lkl-Commander 
and Commander maps. The former results were derived using the 
LklT25693 mask; the latter are those determined previously in 
Sect. 6.4. 




Direction 

Method 

A 

(l,b) [°] 

Lkl-Comm. . , 

. 0.063 ±0.021 

(234, -27) ± 31 

Commander . , 

. 0.067 ±0.023 

(230,-18) ±31 
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7.4. Summary 

Using a larger sky fraction in our analyses leads to small 
changes in the results related to large-angular-scale anoma¬ 
lies, but these are essentially consistent with expectations 
from random Gaussian statistics. In particular, the asym¬ 
metry in power on the sky, as parameterized by a dipole 
modulation model, is robust to mask changes. 

8. Polarization analysis 

As previously discussed in Sect. 2, large angular-scale CMB 
fluctuations in the Planck polarization data have been sup¬ 
pressed by a post-processing high-pass filter to minimize 
the impact of systematic artefacts. Therefore, no polariza¬ 
tion results concerning CMB statistical anomalies on such 
scales are presented in this paper. In addition, a noise mis¬ 
match between simulations and data also limits our ability 
to study polarization more generally. Nevertheless, a local 
analysis of the polarization data for stacked patches of the 
sky can still be performed, in order to test the statistical 
properties of the CMB anisotropies. In this case, the stack¬ 
ing procedure mitigates the impact of the small-scale noise 
and potential systematic effects. 

Traditionally, the Stokes parameters Q and U are used 
to describe the CMB polarization anisotropies (e.g., Zaldar- 
riaga & Seljak 1997). Such quantities are not rotationally 
invariant, thus for the stacking analysis it is convenient to 
consider a local rotation of the Stokes parameters, result¬ 
ing in quantities denoted by Q r and J7 r , as described in 
Sect. 8.1. Additionally, several other related quantities can 
be defined. 

The polarization amplitude P = \J Q 2 + U 2 and po¬ 
larization angle T = ^ arctan (U/Q), are commonly used 
quantities in, for example, Galactic astrophysics. However, 
unbiased estimators of these quantities in the presence 
of anisotropic and/or correlated noise are hard to define 
(Plaszczynski et al. 2014). Of course, a direct comparison 
of the observed (noise-biased) quantity to simulations anal¬ 
ysed in the same manner is possible, but we elect here to 
defer the study of this representation of the polarization 
signal, using maps of the polarization amplitude only to 
define peaks around which stacking can be applied. 

The rotationally invariant quantities referred to as E 
and B modes are commonly used for the global analysis of 
CMB data. Although the E-mode maps are not analysed in 
detail here, they are considered qualitatively, so that it is 
appropriate to recall their construction. Since the quantities 
Q±iU , defined relative to the direction vectors n, transform 
as spin-2 variables under rotations around the ft axis, they 
can be expanded as 

oo i 

(Q ± iU){n) = Y, E ±2 Y tm (h), (62) 

1=2 m=-l 

where ± 2 Wm(^) are the spin-weighted spherical harmonics 

and are the corresponding harmonic coefficients. If we 
define 

a Fm = \ (4m + 4m } ) - ( 63 ) 

a tm = Y (4m - 4m } ) . ( 64 ) 


then the invariant quantities are given by 


00 i 


E(n) 

1=2 m=-l 

00 i 

(65) 

B(h) 

— ^ ^ ^ ^ a £mYlrri(n) • 

1=2 m=-i 

(66) 


8.1. Stacking around temperature hot and cold spots 

The stacking of CMB anisotropies around peaks (hot and 
cold spots) on the sky yields characteristic temperature 
and polarization patterns that contain valuable information 
about the physics of recombination (Komatsu et al. 2011). 
Statistical analysis of stacked images differs from the other 
tests in this paper in several respects. First, peak-related 
new physics may be revealed that is difficult to find in a 
global analysis, for example, the non-Gaussian CMB cold 
spots predicted by a modulated preheating model (Bond 
et al. 2009). Secondly, stacking is a local operation, which 
naturally avoids mask-induced complications. Thus stack¬ 
ing can be used as a transparent and intuitive method to 
test the robustness of anomalies found with other methods. 
Alternatively, it can be applied as a quality indicator of the 
data at the map level. 

Our stacking procedure is as follows. Hot (or cold) peaks 
are selected in the temperature map as local extrema with 
negative (or positive) second derivatives, and classified rel¬ 
ative to a given threshold v (in rms units of the temper¬ 
ature map). Since the spinorial components Q and U are 
expressed in a local coordinate system, we employ a config¬ 
uration in which the Stokes parameters around a peak at 
the direction no can be superposed (Kamionkowski et al. 
1997). In particular, we use a locally defined rotation of the 
Stokes parameters that is written as: 

Q r (n; no) = — Q (n) cos (20) — U (ft) sin (2</>), (67) 

U Y (n; no) = Q (ft) sin (20) — U (n) cos (2</>), (68) 

where p is the angle between the axis aligned along a merid¬ 
ian (pointing to the south by convention) in the local coor¬ 
dinate system centred on a peak at no and the great circle 
connecting this peak to a position n. This definition de¬ 
composes the linear polarization into radial (Q r > 0) and 
tangential (Q r < 0) contributions around the peaks. This 
definition of Q r is equivalent to the “tangential shear” used 
in weak lensing studies. 

For visualization purposes, a fiat patch around each 
peak is then extracted, and the average stacked image com¬ 
puted from the subset. A position on the sky at an angular 
distance 0 from the central peak is labelled with the fiat-sky 
coordinates 

x = w cos <p , y = w sin <p . (69) 

Here w = 2sin(#/2) « 0 is the effective flat-sky radius. For 
the angular scales of a few degrees considered in the stack¬ 
ing analyses the difference between w and 6 is negligible. 
We use w for analyses in the flat-sky approximation, and 0 
for analyses directly on the sphere. 

The stacking process tends to provide an image with 
azimuthal symmetry about its centre, due to the almost 
uncorrelated orientations of the temperature peaks. The 
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stacked images of temperature patches around hot spots se¬ 
lected above the null threshold for both the Commander data 
and a corresponding simulation are shown in the top row of 
Fig. 41. The observed patterns are in excellent agreement. 
Stacking around cold spots yields similar patterns but with 
flipped sign. Given the symmetry, it is often useful to con¬ 
sider the radial profile obtained by averaging the stacked 
image over the azimuthal angle 0. Fig. 42 shows such a 
profile determined from the stacked temperature image. 

At this point, it is useful to consider the underlying 
physics represented by the various patterns in the stacked 
images. During recombination, the sound horizon extends 
an angle 0 S = v s /Da ~ 0.011 (0.61°), where r s « 0.15 Gpc 
is the size of the sound horizon at recombination and 
Da ~ 14 Gpc is the angular-diameter distance to the last 
scattering surface. To understand the ring patterns in the 
stacking image, projection effects must be taken into ac¬ 
count. Firstly, all 3D modes with wavenumber k > if Da 
contribute to a 2D Gmode. More modes contribute to, 
and therefore enhance the power at lower t. For the first 
acoustic peak, the net effect is a 7r/4 phase shift towards 
lower i, such that £ s ~ (tt — 7r/4)/# s ~ 220. The pro¬ 
jected acoustic scale on the temperature map is of order 
# 2D = 7r/4 = 0.014 (0.81°). Secondly, the stacked 2D 
modes around peaks interfere with each other. The first 
dark ring appears at 1.22# 2D « 0.017 (1.0°). The factor 
1.22 is the ratio of the first minimum of the projection ker¬ 
nel, the Bessel function Jo, to the first minimum of the 
unprojected cosine wave. 

The dark ring can also be regarded as a consequence 
of the correlation between T and — V 2 T. At the tempera¬ 
ture maxima —V 2 T is positive, with an amplitude of or¬ 
der T pea k/(# 2D ) 2 . Thus, the quadratic terms in the local 
expansion of the temperature field have a negative contri¬ 
bution that grows as — T pea ^(Ti7 / 0 2B ) 2 . At w > 0 2B the 
quadratic terms dominate and the T-(—V 2 T) correlation 
becomes negative. Meanwhile, the T-(—V 2 T) correlation 
tends to zero on the scale w > # 2D , where the temperature 
autocorrelation becomes weak and the local quadratic ex¬ 
pansion starts to fail. As shown in Fig. 42, the dark ring 
appears at the critical point where the T-(—V 2 T) correla¬ 
tion reaches its minimum and turns back toward zero. 

We have discussed the projection effects that make the 
projected radial acoustic scale on a stacked T image larger 
than 0 S . For Q r , the most striking patterns in the image 
have more intuitive simple explanations, since the stacking 
is essentially the real-space equivalent of the temperature 
polarization correlation. The projection function contains 
an extra f 2 factor, which enhances the high-£ power and 
reduces the projected radial acoustic scale, coincidentally, 
back to « 6 S . The quadrupole responsible for the polariza¬ 
tion around peaks is induced by gravity on angular scales 
larger than twice the size of the horizon at decoupling. In 
the case of an overdensity, this causes a flow of photons 
towards the gravitational well on these scales, inducing a 
quadrupolar pattern (see, e.g. Coulson et al. 1994). The 
spherical symmetry of the gravitational interaction causes 
a rotation of the quadrupole in the vicinity of the well, re¬ 
sulting in a radial configuration in polarization. This radial 
polarization pattern implies Q r > 0 and an overdensity im¬ 
plies T < 0 by the Sachs-Wolfe formulae, which leads to 
anticorrelation on these scales. Similarily, an underdensity 
leads to an outward flow and induces a tangential polariza¬ 
tion pattern, once again leading to anticorrelation on these 
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Fig. 41. From top to bottom, T, Q, U, Q r , and U r stacked 
images (in pK units) extracted around temperature hot spots 
selected above the null threshold (v — 0) in the Commander sky 
map for data (left column) and an equivalent simulation (right 
column). The horizontal and vertical axes of the flat-sky projec¬ 
tion are labelled in degrees. 
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Fig. 42. Radial profile /jLt(vc j) derived from the stacked temper¬ 
ature image (see Fig. 41 or 45). The denominators <jq and <72 are 
the theoretical rms values of CMB T and V 2 T, respectively. The 
theoretical (ijlt(vj)) is a linear combination of (T(tx7)(T(0)/cro)) 
(green dash-dotted line) and {T(vj){— V 2 T(0))/cr2)) (blue dot¬ 
ted line). For all four component-separated maps, the deviation 
of /it from the ensemble mean (/it) of the fiducial model (here 
the Planck 2015 ACDM best fit) is consistent with cosmic vari¬ 
ance, and can be related to the low-^ power deficit. The example 
power-deficit {/a t) (purple dashed line) is the theoretical predic¬ 
tion of (/ it) if the fiducial model Cg s are reduced by 10 % in the 
range 2 < £ < 50. 


scales. At smaller scales, the polarized contribution is domi¬ 
nated by the dynamics of the photon fluid. The acoustic os¬ 
cillations modulate the polarization pattern, leading to the 
different rings in the stacked images. The most noticeable 
rings in the stacked Q r image are approximately at 0 S and 
20 s . Thanks to the £ 2 enhancement, multiple acoustic peaks 
in the TE power spectrum may be captured and projected 
into ring patterns in the stacked polarization images. As 
photons flow towards the overdensity, they are compressed 
and the temperature increases, slowing the fluid descent 
into the well. Eventually, the radiation pressure becomes 
large enough to reverse the photon flow. This expansion 
cools the photons until they fall back towards the well. Note 
that the resulting inner ring was not observed in the WMAP 
analysis (Komatsu et al. 2011), since the resolution was too 
low. 

Figure 41 clearly reveals all of the features described 
above. The two bright rings at 9 S « 0.011 (0.6°) and 
2 0 S « 0.021 (1.2°) are the predicted patterns associated 
with the first Cj B acoustic peak at £ « 310, while the 
two faint rings are a striking illustration of the detection 
of multiple acoustic peaks in the TE power spectrum. The 
large-scale anticorrelation is suppressed due to the scale- 
dependent bias which results from the fact that peaks are 
defined by the second derivatives of the temperature field 
(e.g., Desjacques 2008). 

We are now in a position to discuss the consistency of 
the Planck results with the predictions of a ACDM cosmol¬ 
ogy. For simplicity, further analysis is focused on the angu¬ 
lar profiles, and specifically the mean, /x(0), estimated as the 
average of the angular profiles around all hot (cold) peaks 
above (below) a certain threshold v. This analysis is per¬ 


formed directly on the sphere to avoid any repixelization er¬ 
ror. Note that the expected value of the mean temperature 
angular profile is proportional to f £d£C BT Jo(£@), whilst 
the expected values of the Q r and U v mean angular pro¬ 
files are approximately proportional to f — £d£ Cj E J2{£0) 
and f —£d£Cj B J 2 (£Q\ respectively. Since T has even par¬ 
ity and B has odd parity, the expectation value for Cf B is 
zero, and the U r mean angular profile is therefore expected 
to vanish. 

A x 2 estimator is used to quantify the differences be¬ 
tween the profiles obtained from the data and the expected 
values estimated with simulations: 

x 2 = w) - m\ c- 1 w) - mV , (to) 

with the covariance matrix defined as 

1 N 

E ~ M#*)] K(%) - M(0j)]> ( 71 ) 

1 k =1 

where the sum is over the N simulations used to estimate 
this matrix and jl(0) is the ensemble average. Note that 
although the profiles in Fig. 41 are derived from data at a 
resolution 7V S id e = 1024, faster convergence of the x 2 statis¬ 
tic is achieved using maps at a lower resolution. We have 
verified that the results remain unchanged when adopting 
data with N s [^ e = 512. 

Figure 43 presents a comparison between the profiles ob¬ 
tained from the component-separated data and the mean 
value estimated from simulations processed through the 
SEVEM pipeline. Note that the error bars for the temper¬ 
ature profiles are asymmetric due to a bias in the selection 
of the peaks above a given threshold. Results for hot and 
cold spots are shown for two different thresholds, v = 0 
and v = 3. There is generally excellent agreement between 
the different component-separation methods. A systematic 
deviation between the data and the simulations for the hot 
peaks in temperature {y = 0) is seen at a level greater than 
1 a. This discrepancy increases at higher thresholds, reach¬ 
ing values of about 2 a for the v = 3 case. Similar behaviour 
is seen for the cold spots. For the Q r angular profiles, the 
most striking differences appear around 0 = 2° in the v = 3 
case for hot peaks, and around 6 = 1?5 for the cold peaks. 
For the U r angular profiles, where a null signal is expected 
(i.e., only noise is expected to be present), deviations at 
similar levels are seen. 

Table 38 presents the corresponding p-values for this 
comparison. A theoretical x 2 distribution is used to deter¬ 
mine the probability that a sky drawn from the ACDM 
cosmology has a value larger than that derived from the 
data. We have verified this approach by comparing the em¬ 
pirical x 2 distribution estimated from 100 simulations (in 
which the mean value and the covariance matrix are com¬ 
puted from a further 900 simulations) with the theoretical 
distribution with the corresponding degrees of freedom (see 
Fig. 44). The x 2 value of the data is then estimated using 
the mean value and the covariance matrix determined from 
simulations. Although some differences are found among 
the component-separation methods, a general consistency 
between model and data is found. 

Although the x 2 test has the advantage of being sen¬ 
sitive to different types of deviations between model and 
data, does not assume prior knowledge about possible de¬ 
partures from the model, and can account for correlations 
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Fig. 43. Mean radial profiles of T, Q r , and U T in gK obtained for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). 
Each individual panel contains (top) the mean radial profiles and (bottom) the differences (denoted “Diff”) between the mean 
profiles of the data and those computed from the ensemble mean of the simulations. Results based on stacks around temperature 
hot and cold spots are shown in the left and right columns, respectively. Upper plots present results for peaks selected above the 
null threshold, while lower plots show the equivalent results for peak amplitudes above (hot spots) or below (cold spots) 3 times 
the dispersion of the temperature map. The black dots (connected by dashed lines) depict the mean profiles and the shaded regions 
correspond to the la (68%) and 2 a (95%) error bars. The mean profiles and error bars are determined from SEVEM simulations. 
Note that the Diff curves for each component-separation method are computed using the corresponding ensemble average, although 
only the ensemble average from SEVEM is shown here. 


between the various tests from which it is constructed, it 
can nevertheless be suboptimal under certain conditions. 
This appears to be the case when considering the system¬ 
atic shift between data and simulations seen in the temper¬ 
ature profiles pr — the x 2 statistic is not particularly sensi¬ 
tive to systematic deviations of constant sign. We therefore 
consider an alternative quantity, the integrated profile de¬ 
viation, defined as 

a mw) = [ [fjL T (0) - M0)\ w( 0 ) ae , (72) 

Jo 

where i?, the size of stacking patches, is taken to be 3?5 
in this case. The weighting function is chosen to be pro¬ 
portional to the expected profile, but the results are robust 
for other choices, e.g., W = 1. The p -values obtained in 
this case are given in Table 39. These are consistent with 
what might be expected from visual inspection of the plots, 
i.e., the deviations are typically close to 2 a. These devia¬ 
tions are likely to be connected to the deficit in the ob¬ 
served power spectrum at low multipoles, as may be seen 
in Fig. 42. Here, the purple dashed line indicates the reduc¬ 
tion in JIt if the theoretical Cg values are reduced by 10 % 
over the range 2 < f < 50. 

8.2. Generalized stacking 

In this section, a much wider class of stacking methods is 
introduced, with particular emphasis on oriented stacking , 
a novel approach that has not previoulsy been explored in 
the literature. We regard the stacking as oriented if the ori¬ 
entation of the local coordinate frame, and in particular the 



Fig. 44. x 2 distributions obtained from the T (left column), Q r 
(middle column), and U v (right column) mean radial profiles cen¬ 
tred on temperature hot spots selected above the null threshold 
(upper row) and three times the dispersion of the map (bottom 
row). The black lines correspond to the theoretical % 2 distribu¬ 
tion with 12 degrees of freedom, whilst the histograms show 
the distributions determined from 100 simulations computed 
through the Commander (red), NILC (orange), SEVEM (green), and 
SMICA (blue) pipelines. The vertical lines represent the % 2 values 
obtained from the data. 

(j) = 0 axis, is correlated with the map that is being stacked. 
Thus, the stacking methodology in Sect. 8.1 is considered 
unoriented, because the orientation is defined relative to the 
local meridian pointing towards the Galactic south, rather 
than any property of the data themselves. Alternative ap¬ 
proaches to unoriented stacking can also be considered. In 
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Table 38. p -values of the T, Q r , and U r angular profiles com¬ 
puted from the stacking of hot and cold spots selected above the 
v — 0 and v — 3 thresholds. 




Probability [%] 


Comm. 

NILC 

SEVEM 

SMICA 


v — 0 (hot 

; spots) 



T. 

. 8 

3 

7 

5 

Qr .... 

. 4 

5 

3 

3 

Ur .... 

. 93 

28 

75 

44 


v — 3 (hot 

; spots) 



T . 

. 18 

16 

22 

21 

Qr • • ■ ■ 

. 34 

23 

31 

19 

Ur .... 

. 28 

61 

21 

50 


v — 0 (cold spots) 



T . 

. 23 

38 

29 

39 

Qr .... 

. 86 

85 

63 

78 

Ur .... 

. 14 

11 

39 

34 


v — 3 (cold spots) 



T . 

. 24 

21 

23 

20 

Qr .... 

. 21 

51 

29 

52 

Ur .... 

. 30 

13 

30 

8 


Table 39. p -values of A/it computed from the stacking of hot 
and cold spots selected above the v — 0 and v — 3 thresholds 
from the Commander, NILC, SEVEM, and SMICA maps. 


Probability [%\ 

Comm. NILC SEVEM SMICA 


Hot spots 

T 0 = 0) . 96.0 95.8 96.2 97.1 

T(v = 3) . 98.6 98.2 98.3 98.7 

Cold spots 

T (u = 0) . 97.1 96.9 98.1 97.9 

T(v = 3) . 92.0 90.6 90.6 93.0 


this subsection, the orientation of each patch is chosen ran¬ 
domly from a uniform distribution in [0,27r). The unori¬ 
ented T and Q r images can then be directly compared with 
previous sections. 

For unoriented stacking, the ensemble average of stacked 
fields cannot result in any intrinsic ^-dependence, as this 
would be averaged out by the uncorrelated orientation 
choices. The ^-dependence due to a specific choice of rep¬ 
resentation can always be removed via a local rotation. For 
example, the ensemble averages of Q + iU around unori¬ 
ented temperature peaks are proportional to e 2 ^. A local 
rotation (Q,U) (Q r ,U r ) (Kamionkowski et al. 1997) re¬ 
moves the e 2 ^ factor and compresses the information into a 
single real map Q r . For oriented stacking, the 0-dependence 
can be a mixture of a few Fourier modes (e 2m< ^, for integer 
m). Each m mode corresponds to a radial (^-dependent) 
function. 

In what follows, we use the N s ^ e = 1024 component- 
separated maps at a resolution of 10' FWHM. The use of 
this higher resolution as compared to the N si ^ e = 512 data 
used in Sect. 8.1 is motivated by the smaller-scale features 
that are expected to result from the oriented stacking. 

We also introduce the concept of the noise-free ensemble 
average (NFEA), which is defined as the ensemble average 
of stacked CMB-only maps for a fiducial cosmology. Re¬ 
call that the fiducial model for the simulated sky maps, the 
Planck 2013 best-fit ACDM model (Planck Collaboration 


XVI 2014), differs from the updated Planck 2015 best-fit 
ACDM model (Planck Collaboration XIII 2015). In previ¬ 
ous sections, this mismatch was partially accommodated 
by rescaling the CMB signal by a fixed scale factor. Here, 
we instead specifically adopt the 2015 best fit as a fiducial 
model for the data. When comparing the data to the simu¬ 
lations, we subtract the corresponding NFEA to minimize 
any bias resulting from cosmology dependence. 

In the context of random Gaussian fields, the NFEA 
can be computed straightforwardly following Bardeen et al. 
(1986): 

(M) = ( k Mw t }(ww t ) (73) 

where M is the map (around the central peak) to be 
stacked, and w is the collection of Gaussian variables (on 
the central peak) that are related to peak selection and ori¬ 
entation determination. Eq. (73) is only valid for Gaussian 
random variables. If the patch is rotated before stacking, 
the field value evaluated at a dynamic coordinate is, in gen¬ 
eral, not a random Gaussian variable. However, statistical 
isotropy guarantees that the rotation of patches is equiva¬ 
lent to an orientation constraint on the nonzero-spin field. 
For example, orienting each patch in the direction where 
U = 0 and Q > 0 is equivalent to the unoriented stacking 
case where only peaks satisfying the additional constraint 
—e/2 < arg (Q + iU) < e/2 (e —» 0 + ) are selected. 

A further source of statistical bias can arise from noise 
mismatch between the simulations and the data. Since the 
effect of noise is to introduce random shifts in the peaks and 
hence suppress patterns in the stacked images, any noise 
mismatch can lead to pattern mismatch between the data 
and simulations. For the temperature data, the contribution 
due to noise mismatch is estimated to be at the sub-percent 
level, lower than the cosmic variance. For stacking on polar¬ 
ization peaks, the impact of the noise mismatch cannot be 
safely ignored. Thus, for quantitative comparisons in this 
paper, we only consider stacking on temperature peaks. 


8.2.1. Oriented temperature stacking 

The most straightforward way to orient a patch centred 
on a temperature peak is to align the horizontal axis with 
the major axis defined by a local quadratic expansion of 
the temperature field around the peak. The disadvantage 
of doing so is that the orientation is dominated by small- 
scale fluctuations that are noise-sensitive. A better choice 
is to use the major axis of the inverse Laplacian V -2 T that 
filters out the small-scale power. The inverse Laplacian is 
defined as: 


rt i») = -EE ( wti) y, ~ (ft) • (74) 

where aj m are the harmonic coefficients of the masked tem¬ 
perature map. Spin-2 maps Ut are then defined by: 

oo i 

(Qt ± iU T ) (ft) = £ X «L [ ±2 YUn)} • (75) 

1=2 m=—l 

In the flat-sky limit, Qt ~ (d 2 — <9 2 )(V -2 T) and Ut ~ 
—2 d x d y \7~ 2 T. To align the V -2 T axes of the patches, we 
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rocos <f> 


wcos <fi 


Fig. 45. Comparison between unoriented stacking (upper pan¬ 
els) and oriented stacking (lower panels) of temperature patches 
around temperature hot spots selected above the null threshold 
(v — 0). The left panels are the stacked SMICA maps, and the 
right panels their corresponding NFEAs. The image units are 
fiK. 


rotate each patch so that Ut vanishes and Qt > 0 for the 
central peak. 

Figure 45 presents the stacked images of SMICA temper¬ 
ature patches centred on temperature hot spots selected 
above the threshold v = 0, in both unoriented and oriented 
forms. These are seen to be in excellent agreement with 
their accompanying NFEAs, and, in the case of the unori¬ 
ented stacks, with the results shown in Fig. 41, despite the 
different stacking methodologies adopted (and component 
separation method selected for visualization purposes). 

The oriented T image is notably different from the un¬ 
oriented one. The alignment between the major axis (of 
V -2 T) and the horizontal axis results in an ellipse elon¬ 
gated along the horizontal axis, rather than a central disc. 
Moreover, the quadratic-term contribution is suppressed 
along the horizontal axis where the temperature profile is 
smoother, and enhanced along the vertical axis where the 
temperature profile is sharper. As a consequence, the dark 
ring visible in the upper panel at 1° splits into two cold 
blobs along the vertical axis. 

To proceed with a quantitative analysis, we extract 
Fourier modes T m (zu) from the stacked map T stac k(^,0) 
as follows: 

1 f 2n 

T m (w) = - —— / T stack (m, 4>) cos to0 dcj), (76) 

(1 + d m0 )7T J o 

where S m o is the Kronecker delta function. For odd m, the 
NFEA (T m ) vanishes due to statistical isotropy. For even m, 
a straightforward calculation shows that only To(tu), which 
is equivalent to and X^tu) have nonzero NFEAs. 

As discussed previously in Sect. 8.1, there are some 
shortcomings of the standard y 2 procedure that is generally 
used to assess the consistency of the data with simulations. 


The problem is simplified by studying the statistics of an 
integrated profile deviation: 

Tm(W) = [ [T m (w) - (T m (w))\ W[w) dw , (77) 

Jo 

where X, the size of the stacking patches, is taken to be 
2° in our examples. The purpose of removing the NFEA, 
(T m (tu)), which differs for the data and the simulations, 
is to minimize the impact of the cosmology dependence. A 
natural choice for the filter is (T m (w)) itself with a proper 
normalization: 


W(w) 


(T m (w)) 
irf {T m {w)) 2 dw 


(78) 


For the filter given by Eq. (78), the integrated profile devia¬ 
tion Tm describes the relative deviation from the NFEA. If 
ACDM is the correct model, the deviation is due to cosmic 
variance and noise. The distribution of T m is obtained from 
simulations. 

Table 40 presents a comparison of the T m values de¬ 
rived from the Planck data and the FFP8 simulations. No 
inconsistencies in excess of the 3 a level have been found, 
although tensions around 2 a are seen. 

The m = 0 projection kernel Jo[(£ + l/2)w] peaks at 
low t. Thus To is cosmic-variance sensitive and the apparent 
discrepancy in it could be related to a \ow-£ power deficit. 
An example is shown in Fig. 42 for illustration. To test 
the robustness of this result, we have tried three additional 
filters: a top-hat filter W = 1, a linear filter W = w, and 
a Gaussian filter W = exp(— w 2 /a 2 ) with cr g = 1°. In all 
cases, the power deficit remains at about the 2 a level. 

Although the To deficit is not significant enough to fal¬ 
sify the ACDM model, further investigation of its proper¬ 
ties may still be interesting and help us understand the 
other anomalies discussed in this paper. We consider two 
possibilities. Firstly the amplitude of the To deficit is of or¬ 
der 5-10%, which coincides with the level of hemispherical 
power asymmetry discussed in Sect. 6.1. To test whether 
the To deficit is localized on one hemisphere, we define the 
“north” direction to be aligned with the power asymmetry 
direction at (Z, 6) = (212°,—13°) (Akrami et al. 2014) and 
compute To on the northern and southern hemispheres sep¬ 
arately. The results are presented in Table 41. Although the 
To deficit is more significant for the southern hemisphere, it 
remains consistent with the ACDM prediction. Secondly, it 
is of interest to determine whether the To deficit is related 
to the Cold Spot discussed in Sect. 5.7. We therefore mask 
out the Cold Spot using a disc of radius 6° and repeat the 
calculation. The impact of this region on the To deficit is 
insignificant. 

Tensions at the 2 a level are also seen for %. However, 
due to the additional i 2 factor in the projection kernel, the 
oriented (m = 2) component T is more sensitive to high-£ 
power where the cosmic variance is small, and an under¬ 
standing of the noise properties of the data is more im¬ 
portant. The former implies that the related uncertainty in 
T is, in general, smaller than that in To • However, a mis¬ 
matched cosmology, perhaps arising from a different pri¬ 
mordial power amplitude A s , can then lead to significant 
tension between the data and the simulations. Indeed, we 
find that without application of our cosmology calibration 
(i.e., the subtraction of the NFEA in Eq. 77) the ^-tension 
between the data and simulations increases by about 0.5 cr, 
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Table 40. Tm , as defined in Eqs. (77) and (78), for different thresholds v. The expected values, together with the la (68% CL) 
and 2 a (95 % CL) ranges determined from simulations are given in brackets. 




To 


r 2 

Method 

Hot spots 


Cold spots 

Hot spots Cold spots 


threshold v — 0 


Commander . . 

1 

0 

b 

^0 

rv rj a +0.06+0.15 \ 
U,U ^-0.04-0.07f 

1 

O 

b 

rv +0.06+0.16 \ 

U,U ^-0.04-0.07f 

0.07(0.04!?;??!?;??) 

0.06(0.04!?;??!?;?? 

NILC . 

. . —0.04( 

n c\A +0.06+0.15 \ 
U,U ^-0.04-0.07f 

—0.05( 

n r\A +0.06+0.16 \ 
u,u ^-0.04-0.07f 

0.06(0.03!?;??!?;??) 

0.05(0.03!?;??!?;?? 

SEVEM . 

. . —0.03( 

n c\A +0.06+0.16 \ 
U,U ^-0.04-0.07f 

—0.04( 

n r\A +0.06+0.16 \ 
u,u ^-0.04-0.07f 

0.06(0.04!?;??!?;??) 

0.06(0.04!?;??!?;?? 

SMICA . 

. . —0.03(- 

n m +0-03+0.07 \ 

u,UJ --0.02-0.04J 

—0.05(- 

r» nn+0.03+0.07\ 
U.UU- 0 .o2-0.04f 

0.06(0.04!?;??!?;??) 

0.06(0.04!?;??!?;?? 





threshold v 

= 1 


Commander . . 

. . —0.06( 

r» nir+0.09+0.22 \ 
U - UO -0.05-0.10f 

—0.06( 

n n^^^-oo+o^i ^ 
U-UiJ_o.o6-0.09/ 

0.06(0.03!?;??!?;??) 

0.04(0.03!?;??!?;?? 

NILC . 

. . —0.06( 

r» nir+0.09+0.22 \ 
U - UO -0.05-0.10f 

—0.07( 

n n^^^-oo+o^i \ 
U.UU-o.05-0.09/ 

0.06(0.02!?;??!?;??) 

0.04(0.02!?;??!?;?? 

SEVEM . 

. . —0.06( 

n oc _ i _ o-09+o.22 \ 
U - UO -0.05-0.10f 

—0.06( 

n n^H^-09+0.22 \ 
U - UO -0.05-0.10J 

0.06(0.03!?;??!?;??) 

0.04(0.03!?;??!?;?? 

SMICA . 

. . —0.06(- 

n ni +0.04+0.10 n 
u - u± -0.03-0.06/ 

—0.07(- 

n ni +0.04+0.10\ 
UlUJ --0.03-0.06/ 

0.06(0.03!?;??!?;??) 

0.04(0.03!?;??!?;?? 


Table 41. 7o, as defined in Eqs. (77) and (78), for different thresholds v and hemispheres. The “north” hemisphere is centred on 
the Galactic coordinate (l, b) = (212°, —13°) and the “south” hemisphere in the opposite direction. The expected values, together 
with the 1 a (68 % CL) and 2 a (95 % CL) ranges determined from simulations are given in brackets. 


Method 


Commander 

NILC 

SEVEM . . . 
SMICA . . . 


Commander 

NILC 

SEVEM . . . 
SMICA . . . 


“North” To 


“South” 7o 


Hot spots 


—0.02( 0.03!?;??!?;??) 
— 0 - 02 ( 0 . 021 °;™?) 
-0.02( 0.03!?;??!?;??) 
- 0 . 02 + 0 . 01 !?;??!?;??) 


—0.04( 0.03!?;??!?+) 
—0-05( 0.03!?;??!?+) 
—0.04( 0.04!?;??!?+) 
-0.04(-0.02!?;??!?;? 3 7 ) 


Cold spots 


Hot spots 


-0.03( 

-0.03( 

-0.03( 


threshold v — 0 


n nq+°- 07 +°- 18> i 
U.UO_ 00 4-0.07f 

n n9 +0.07+0.17\ 
u.uz, —0 04 —0.07/ 


J —0.04 —0.07 
i +0.04+0.09 


—0.05( 
—0.05( 
—0.05( 


n nQ+°- 07 +°- 18> \ 
U.UO-o 05-0.07f 

n n9 +0.07+0.18\ 
U.UZ_o 04-0.07/ 


-0.03(-0.0l!?;?™) —0.05(—0.0l!?;??!?;??) 


-0.05-0.07 

+0.04+0.08 


threshold v — 1 


—o.o5( o.o3±S:^±S:?S) 
—0-06( 0.02+;??!?+) 
—0-05( 0.03+;??!?+) 
-0.05(-0.02+?t+;? 3 7 ) 


—0-08( 0.041?;™?) 
—0-08( 0.031?;??!?;??) 
—0-08( 0.04!?;??!?;??) 
-0.08+0.02+;??+;??) 


Cold spots 


-0.06( 0.03+;??+;??) 
-0.06( 0.02!?;?™) 
—0-06( 0.03+;??+;??) 
-0.07+0.01+;??+;??) 


—0-08( 0.04!?;??!?;??) 
—0.08( 0.031?;??!?;??) 
—0-08( 0.04!?;??!?;??) 
—0.09(—0.02!?;??!?;??) 


whereas the variation of the 7o-tension is < 0.2 a. The high- 
£ sensitivity of T 2 also requires the use of an accurate noise 
model, and it is possible that the 1-2 a tension in T 2 may 
be alleviated once improved noise simulations are available. 

8.2.2. Oriented polarization stacking 

The stacked Q and U images can be decomposed into 
Fourier modes, Q + iU = Ylm=-oo For unori¬ 

ented Q + iU stacking on temperature peaks, only P 2 (tu) 
has a non-zero NFEA, and it can be linked to the con¬ 
ventional Q r stacking via P 2 = — Q r . Figure 46 shows 
that the stacked Q r image is in excellent agreement with 
its NFEA and the corresponding stacked image (fourth 
panel) in Fig. 41, despite the different stacking methodolo¬ 
gies adopted (and component-separation method selected 
for visualization purposes). The length and orientation of 
the headless vectors represent the polarization amplitude, 

p stack = VQLck + Gtack> and direction. 

We next consider oriented stacking of the polarization 
maps, again using Qt, Pt to define the orientation of the 
patches. The stacked polarization images around temper¬ 
ature peaks have m = 0,2,4 Fourier components. We can 



Fig. 46. Stacked Q r image around temperature hot spots se¬ 
lected above the null threshold (v — 0) in the SMICA sky map. 
The left panel corresponds to the observed data and the right 
panel shows the NFEA. The image units are /iK. The head¬ 
less vectors (black solid lines) are the polarization directions for 
stacked Qstack, Fstack- The lengths of the headless vectors are 
proportional to the polarization amplitude P s tack- 


also choose to stack the polarization maps on Pt peaks, 
where Pt = \JQ\ + P|>. This picks up m = 0,4 Fourier 
modes with no circularly symmetric (Q r , m = 2) mode. In 
Fig. 47 we compare the (Q, U) images stacked centred either 


Article number, page 55 of 66 






















































A&A proofs: manuscript no. planck_2015_iands_final 



I ' I • ! ' I ! I 1 ! « I | l ' I 1 I • ! ' I ' l ' r 



_ I --j. , . . - .4.._.. . L - j I .._ i _ 1 

-0.025 0 0.025 -0.025 0 0.025 

VJ COS (f) ZU COS (j) 




I 


I 

I 

i 


LO 

o 


o 


LO 

o 

I 


o 


Fig. 47. Oriented stacking of polarization fields (Q, U) on tem¬ 
perature maxima (upper panels) and Pt maxima (lower panels). 
In both cases the threshold v — 0 is used and the orientation is 
chosen such that Ut — 0 and Qt > 0 on the central peak. The 
image units are pK. The left panels are the stacked SMICA maps, 
and the right panels their NFEAs. See Fig. 46 for the meaning 
of the headless vectors (black dashed lines). 


Fig. 48. Top : E-mode maps stacked on the unoriented E-mode 
maxima computed above the null threshold v — 0. Bottom: Q 
stacked around oriented polarization amplitude (P) maxima. In 
this case, no threshold is used and the orientation is chosen such 
that U — 0 and Q > 0 on the central peak. The left panels are 
the stacked SMICA maps, and the right panels their corresponding 
NFEAs. See Fig. 46 for the meaning of the headless vectors 
(black dashed lines). The image units are pK. 


on T peaks (top panel) or on P T peaks (bottom panel) with 
their corresponding NFEAs, and find excellent agreement. 

For a quantitative comparison, we only consider stack¬ 
ing on temperature peaks and define the polarization inte¬ 
grated profile deviation 


r R 

V m (W) = / (P m [w) ~ <P ro (07)» W{w) dw , 
Jo 


where by default the filter is 


W{w) 


( P m { vj)) 

/<f (PnM ) 2 dm 


(79) 


(80) 


The comparison of P m (m = 0,2,4) between the data and 
the simulations is shown in Table 42, where the results are 
seen to be in excellent agreement. 

Finally, we note that the peak selection does not have 
to be made from the temperature map. In Fig. 48 we show 
a few examples of stacking on polarization peaks using the 
TVgide = 512 maps. The higher-resolution polarization data 
are too noisy for peak selection. In the upper panels, we 
compare stacked images of the F-mode map centred around 
E-mode peaks with the corresponding NFEA. We find that 
the noise impact is relatively minor for FWHM = 20' maps 
and the plots are in qualitatively good agreement. Another 
possibility, shown in the lower panels, is to stack polar¬ 
ization maps centred on peaks determined from the cor¬ 
responding polarization amplitude map. In this case the 
peaks are strongly biased by the quadratic noise contribu¬ 
tion and quite visible deviation from the NFEA is observed 
in the stacked image. 


9. Conclusions 

In this paper, we have presented a study of the statisti¬ 
cal isotropy and Gaussianity of the CMB using the Planck 
2015 data, including the full mission for temperature. We 
do not claim that our results support or refute any partic¬ 
ular physical model. Rather, we focus on null-hypothesis 
testing: a number of tests are performed, then p -values are 
calculated and reported. It is in the very nature of such a 
model-independent approach to leave the detailed interpre¬ 
tation to the reader. However, we do address the important 
subject of a posteriori correction where possible. 

The statistical tests are performed on maps of the CMB 
anisotropy that result from the application of the four 
component-separation methods described in Planck Col¬ 
laboration IX (2015). All of the results presented here are 
robust with respect to the choice of component-separated 
CMB map. This is important since it demonstrates the high 
quality and equivalence of the Planck component-separated 
data products rendered by different methodologies under 
varying assumptions. 

We find that the CMB is largely consistent with statisti¬ 
cal isotropy, although there are a few indications of anoma¬ 
lies with respect to the expectations of ACDM. Some of the 
tests we have performed are the same as those in PCIS13, in 
which case the results are consistent. Since many of these 
anomalies were also observed in the WMAP temperature 
data, we re-emphasize explicitly the statement we made in 
2013 — that the agreement between the two independent 
experiments effectively rules out the possibility that the 
origin of these features can be found in residual systematic 
artefacts present in either data set. We have also performed 
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Table 42. Pm-, as defined in Eqs. (79) and (80), for different thresholds v. The expected values, together with 1 cr (68% CL) and 
2 a (95 % CL) ranges determined from simulations are given in brackets. 


Method 

Vo 

V 2 

Va 


Hot spots, threshold v — 

0 

Commander . . . 

• 0.06(0.02!° 0 ;“!° 0 ° 5 6 ) 

—0.01( O.Ol!™™) 

0.04(0.01!°;°™) 

NILC . 

. 0.05(0.02+H3 + o.o6) 

—0.02( O.OOl™™) 

0.03(0.0ll° 0 ;° o ™) 

SEVEM . 

. 0.05(0.02 +0 o ;° 3 3+° o ° 6 6 ) 

0 . 01 ( o.oi!™™) 

0.04(0.01!™™) 

SMICA . 

. 0.05(0.03 + °;° 3 3 + o.o5) 

-0.02( 0.001°;°™) 

0.03(0.01!™™) 


Cold spots, threshold v — 

0 

Commander . . . 

. 0.06(0.02i“;“t°;“) 

-0.01( 0.01!°;°™) 

0.03(0.01!™™) 

NILC . 

. 0.06(0.02l°;° 3 l° o ;° 5 6 ) 

-0.01( O.OOl™™) 

0.04(0.011™™) 

SEVEM . 

. 0.06(0.03lS;° 3 3 l o o ;° 5 6 ) 

0.01( 0.0ll°;°il°;° 2 2 ) 

0.03(0.011™™) 

SMICA . 

. 0.05(0.03!™™) 

—0.02( 0.00!°;°il°;° 2 ) 

0.03(0.02!™™) 


Hot spots, threshold v — 

1 

Commander . . . 

. 0.04(0.02l°;° 3 l° o ;“) 

-0.02(-0.00l™™) 

0.05(0.011™™) 

NILC . 

. 0.06(0.02l°;° 3 l° o ° 7 7 ) 

-0.02(-0.0ll°;™ 3 3 ) 

0.05(0.011™™) 

SEVEM . 

. 0.05(0.02l° o ° 4 4 l°;° o 7 ) 

-0.01(-0.00!°;° 2 1°;° o 3 3 ) 

0.05(0.011™™) 

SMICA . 

. 0.04(0.031™™) 

—0.02( 0.00l°;°il°;° 2 ) 

0.06(0.011°;°™) 


Cold spots, threshold v — 

: 1 

Commander . . . 

. 0.07(0.02l°;° 3 l° o ° 6 7 ) 

-0.00(-0.0ll°;° o 2 l°;° 3 ) 

0.01(0.011°;°™) 

NILC . 

. 0.08(0.02l° o ° 3 l™ 6 7 ) 

-0.01(-0.0l!°;™ 3 ) 

0.01(0.011°;°™) 

SEVEM . 

. 0.09(0.02l™ 3 3 l™ 7 ) 

-0.00(-0.00l°;° 2 l°;° 3 3 ) 

0.02(0.011°;°™) 

SMICA . 

. 0.06(0.03 + q'q 3+ qq 7 ) 

—0.01( 0.00!°;°™) 
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a number of new tests, in order to try to narrow down the 
nature of the apparent violations of statistical isotropy. In 
addition, although the component-separated polarization 
maps contained in the Planck 2015 release are high-pass 
filtered, we have performed a stacking analysis that tests 
some aspects of the polarized sky while mitigating the im¬ 
pacts of noise and systematic effects. 

In Sect. 4, we examined aspects of the Gaussianity of 
the CMB fluctuations. Tests of skewness, kurtosis, multi¬ 
normality, 7V-point functions, and Minkowski functionals 
yielded no indications of significant departures from Gaus¬ 
sianity, while the variance of the CMB map was found 
to be low, in agreement with previous studies (PCIS13). 
First-order moments of filtered maps also exhibit the low- 
variance anomaly, as well as a kurtosis excess on certain 
scales associated with the Cold Spot. A new study of peak 
statistics finds results consistent with the expectations for 
a Gaussian random field, although the Cold Spot is again 
detected. 

Section 5 provides an updated study of several previ¬ 
ously known peculiarities. We study in detail the low vari¬ 
ance anomaly, which appears to be associated with the 
known low-£ deficit in the angular power spectrum. We 
confirm the lack of large-scale angular correlations, rela¬ 
tively featureless northern ecliptic hemisphere 3- and 4- 
point functions, and indications of violations of point- and 
mirror-parity symmetry, although we make little or no at¬ 
tempt to correct these for a posteriori effects. We place 
tight constraints on a quadrupolar power modulation. The 
Cold Spot is examined further, and, while we find variance, 
skewness, and kurtosis angular profiles consistent with the 
expectations of statistically isotropic simulations, the mean 
temperature profile is anomalous at roughly the 1 % level, 


apparently due to the surrounding hot ring — the feature 
that deviates most from the Gaussian model. 

In Sect. 6 we perform a series of tests probing the well- 
known large-scale dipolar power asymmetry. We detect the 
asymmetry via pixel-to-pixel variance, as well as by mea¬ 
suring power explicitly or indirectly via i to i db 1 mode 
coupling. The latter approach lends itself to a posteriori 
correction, which reduces the significance of the asymmetry 
substantially when no model for the anomaly is assumed. 
In addition, we perform two independent but related tests 
of directionality. One finds suggestions of anomalous clus¬ 
tering of directions out to relatively small scales while the 
other does not, evidently due to being optimized for slightly 
different forms of directionality. 

Section 7 demonstrates that the significances of several 
large-angular-scale anomalies are robust to the use of larger 
sky coverage, with the observed small changes being con¬ 
sistent with expectations from random Gaussian statistics. 

Finally, Sect. 8 presents the results of the stacking of 
temperature and polarization peaks. We find results that 
are largely consistent with statistically isotropic simula¬ 
tions, both for oriented and unoriented stacking. The excep¬ 
tion is a low unoriented temperature profile, which seems 
to be yet another reflection of the large-scale power deficit. 

With the Planck 2015 release, we are probably near the 
limit of our ability to probe the CMB anomalies with tem¬ 
perature fluctuations alone. The use of large-angular-scale 
polarization, expected for the final Planck release, should 
enable independent tests of these peculiar features. Impor¬ 
tantly, this will reduce or eliminate the subjectivity and am¬ 
biguity in interpreting their statistical significance. It is a 
tantalizing possibility that some of the anomalies described 
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in this paper will take us beyond the standard model of 
cosmology. 
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Appendix A: Generalized Savitzky-Golay 
polynomials 

In the construction of optimal linear filters, one needs 
to combine information about the (statistically isotropic) 
CMB signal, anisotropic instrumental noise, masking to be 
applied for the elimination of foreground contributions, and 
a model for any non-Gaussian signal for matched filtering. 
These can be combined in a general framework of normal¬ 
ized convolutions (Knutsson & Westin 1993), where the fil¬ 
tered field is defined as 


aB ★ wT 
aB ★ B^w ’ 


(A.l) 


where B is the (multiscale) filtering beam function, T is 
the temperature, a and w their respective weights, and * 
denotes the usual convolution operation 


{aB * wT}(£) = ^ a{x)B(x) • — x)T {£ — x). (A.2) 

X 

In the absence of a specific model for the non-Gaussian 
signal, the beam functions can be taken to be orthogonal 
polynomials on a disc, weighted by some smoothing func¬ 
tion, while the weights applied to the temperature maps 
are determined by the CMB and noise covariance. 

In a simple approach, the information about the CMB 
signal can be utilized by pre-whitening the map by convolv¬ 
ing it with an isotropic beam function we = CJ 1 ^ 2 derived 
from the isotropic best-fit CMB power spectrum combined 
with a diagonal approximation to the instrumental noise 
covariance. After the component-separated CMB maps are 
pre-whitened, and the corresponding mask is applied to the 
resulting map, the multiscale filtering kernel be is applied 
at various scales. 

In this paper, the maps are pre-whitened with the 2013 
best-fit cosmological parameter CMB spectrum (Planck 
Collaboration XV 2014), co-added to an isotropic noise 
power spectrum derived from the half-mission, half¬ 
difference noise maps appropriate for each component- 
separation method. No adjustment is made either for the 
recalibration of the 2015 data relative to the nominal re¬ 
sults that the cosmological spectrum is derived from, or 


for the mismatch in noise level between the half-mission, 
half-difference and full-mission maps. This implies that the 
filtering is sub-optimal, but the data and simulations are 
treated consistently so there should be no significant impact 
on the results. The resulting pre-whitening beam function 
we for the SMICA temperature map is shown in Fig. A.l. 

The peak detector wavelets are taken to be Savitzky- 
Golay polynomials (Savitzky & Golay 1964), generalized to 
be defined on a disc with a polynomial smoothing weight 
function applied, as shown in Fig. A.l. A generalized spher¬ 
ical Savitzky-Golay kernel of order n and smoothing weight 
k (referred to as SSGnk in the text) is defined by a polyno¬ 
mial function of a radial variable x = sin(0/2)/ sin(# max /2), 


Fn,k (%) 


' n/2 

X a i x2z 

i =0 


(l-x 2 ) 


2 \k 


(A.3) 


which is normalized to have unit mean on a disc and is 
orthogonal to all non-constant polynomials up to order n, 


i 

/■ 


i 

/ 


xF ntk (x) Ax = 1, / x l+1 F ntk (x) dx = 0. 


(A.4) 


These are essentially high-order low-pass filters in harmonic 
space, but have compact support on the sphere. A few rep¬ 
resentative Savitzky-Golay polynomials are compared to a 
Gaussian kernel in Fig. A.l. Combined with pre-whitening, 
the total effect of the filters applied is described by the 
composite beam functions shown in Fig. A.l. 

One should note a slight Uspace bandwidth mismatch 
between differently shaped kernels with the same FWHM 
value in real space, which is clear from the lower left panel of 
Fig. A.l. While not a problem in general, some care should 
be exercised when directly comparing results for different 
shape kernels. In particular, the £ value at which the filter 
kernel coefficient reaches be = 5 max /2 differs by a factor 
of 1.58 between the GAUSS and SSG84 kernels of the same 
FWHM. 


Appendix B: Doppler boosting 

The main effect of our relative motion with respect to the 
CMB rest frame is a dominant contribution to the CMB 
dipole (Ci); this is boosting of the monopole and has been 
detected previously (Kogut et al. 2003; Fixsen et al. 1996; 
Hinshaw et al. 2009). A subtler consequence of our motion is 
the boosting of all other multipoles. In fact, there are really 
two effects at work. The first is a modulation effect which 
increases power by approximately 0.25 % in the direction 
of our motion and decreases it by the same amount in the 
opposite direction. This can equivalently be thought of as 
coupling between the multipoles i and i ± 1. The second 
is an aberration effect which shifts the apparent direction 
in which CMB photons arrive at our detectors toward the 
velocity direction. 

Planck Collaboration XXVII (2014) reported a detec¬ 
tion of this Doppler boosting, and an associated mea¬ 
surement of its velocity signature of 384 ±78 (statistical) 
± 115 (systematic) kms -1 in the known dipole direction, 
(Z, 6) = (264°, 48°). Here, we demonstrate that the Planck 
2015 data release remains in agreement with this result, by 
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Fig. A.l. Generalized Savitzky-Golay polynomials are orthogonal to polynomials up to degree n on a disc, with smoothing weight 
applied. Upper left panel shows a few representative polynomial kernels (SSG21 in red, SSG42 in dark green, SSG84 in blue) and 
Gaussian (in black) as a function of radius (scaled to the same FWHM of 800 x ), lower left shows their harmonic space representation. 
Right column shows pre-whitening kernel for SMICA temperature map on the top (in light blue), and the corresponding composite 
kernels (WHITE*SSG21, etc) on the bottom (in the same colors). 


Table B.l. Significance measures for the /3 estimates for the 
143 x 217 data set. \ 2 is formed from the three modes of /3 using 
the covariance matrix measured from Doppler boosted simula¬ 
tions. 


Estimator 

X 2 

PTE [%} 

h . 

3.28 

7.01 

Al. 

0.21 

64.39 

K . 

0.08 

77.53 

p . 

3.38 

33.70 


considering the angular scales 500 < £ < 2000. However, 
since the simulations employed in the analysis contain the 
effects of Doppler boosting, we report a consistency check 
rather than a detection. 

It is useful to perform a harmonic transform on the pe¬ 
culiar velocity vector, 

pLM = J dhYZ M (h)/3-n, (B.l) 

where only the L = 1 modes are non-zero. Following the 
convention in Planck Collaboration XXVII (2014), we ro¬ 
tate to an orthonormal basis, labelled (3^ (along the ex¬ 
pected velocity direction), f3 x (parallel to the Galactic 
plane), and /3j_ (the remaining vector). 

The peculiar velocity is detected using estimators that 
pick out the off-diagonal components of the CMB covari¬ 


ance matrix 




2^2 /CMB 


= X(-D 
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LM 
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m\ m 2 M 


(2li + l)(2l 2 + 1)(2£ + 1) 


47r 


£i £ 2 L 


Plm- (B.2) 


The weight function W@ v is a sum of the modulation 
(b v W r ) and aberration (W^) effects. We quote results 
based on orthogonalized weight matrices, 


w* = w (t) - w T n (t)T /n TT (b.3) 

FF f = W T - (B.4) 

Due to the clear connection between the velocity estima¬ 
tors and those used for the lensing analysis, we adopt the 
same data (143 GHz and 217 GHz sky maps, with dust fore¬ 
grounds removed using the 857 GHz data as a template) 
and mask as used in Planck Collaboration XV (2015). The 
results are summarized in Table B.l, and show good con¬ 
sistency with previous results. 


Appendix C: Generalized modulation estimator 

Consider a parameter X that the (primary) CMB power 
spectrum is dependent on. Let X have a dipolar dependence 
of the form X{n) = Xq + A Xh • rfi (this could correspond 
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to a gradient in X across our observable volume), where Xq 
is the average value, ft is the direction to the last scattering 
surface, and m is the gradient direction. To linear order in 
XX/X, the measured spherical harmonics coefficients are 
given by 

dn iso 

= afa + - , (C.l) 

M £'m' 


where the af// are the unmodulated statistically isotropic 
modes. The £Jdm£'m' are coupling coefficients given by 

^>£m£'m' ^rn'rn T &£'t^rl^tm) 5 (^.2) 

^>£m£'m' — ^rn'm^fl {$£'1— \^t— 1 ±m—JL ^£'£-\-l-^£^fm) 5 (0.3) 

where 


A dm — 


-^£m — 


' (£ + l) 2 — m 2 
(2£ + l)(2£ + 3 )’ 

I (£ + m + l){£ + m + 2) 

2(2^ + 1)(2^ + 3) 


(C.4) 


(C.5) 


From Eq. (C.l) we can find the covariance matrix to first 
order in the components A Xm- 


C £m£'m' — CdSu'5 mm' + G^E A X M $. 


M 

'm£'m' ) 


(C.6) 


M 


where 5Cu +1 = dC^/dX + dCt+i/dX. To determine the 
best-fit parameters, we proceed by maximizing the CMB 
likelihood function 


Upon switching bases, we find 



£rn 


Cr(ax +1 ),sr(ax+i) - 2 


E 


C/C/+1 


/4 2 


2 

777 ^ 


( 0 . 12 ) 

(C.13) 


We can then assign the standard errors, <7 = \/E _1 . 


Appendix D: Weighted-variance modified shape 
function estimator 

The BipoSH representation characterizes the off-diagonal 
elements in the covariance matrix and is a generalization of 
the angular power spectrum, Ce, 

j/ ^£ 1 £ 2 = y! ( a £im 1 a £ 2 m 2 )C'£ 1 m 1 £2m2' (B-l) 

m\m 2 

In general, it is not possible to analyse the full sky even for 
component-separated maps, due to the presence of resid¬ 
ual contributions from diffuse Galactic emission and point 
sources. However, the application of a mask leads to cou¬ 
pling between the spherical harmonic modes. Hence, the 
correlation function is no longer described only by C(Q) or 
the power spectrum Ci, and other quantities are required 
to completely quantify the statistical field. 

We obtain an analytic expression for the observed Bi¬ 
poSH coefficients after the application of a mask in terms 
of the corresponding coefficients of the unmasked sky, and 
those of the mask itself, 


V^TrjCj 


exp(— d)C 1 d/ 2), 


(C.7) 


where d is the CMB temperature data. Equation (C.7) is 
maximized for the A Xm that satisfy 


dtC~ 


AC 


AAX M 


-C~ 1 d = Tr 
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d C 


AAX M 


(C.8) 


From Eq. (C.6) it is clear that the CMB covariance 
can be decomposed into an isotropic part (Ci) and a 
small anisotropic part proportional to XX m- By inverting 
Eq. (C.6) and using the orthogonality of the we 

can determine the best-fit parameters 


6 £m Cgct +1 ^£m a £m a £+1 rn 
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(C.9) 

(C.10) 


and Al_i = — AWjG, to first order in the anisotropy. 
These estimators are the full-sky, no-noise versions of 
Eqs. (44) and (45). 

Errors can easily be found by expanding the log- 
likelihood about the best-fit parameters. The Fisher matrix 
is defined as 


Fmm' = -Tr 
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D.2) 


where = y/2£ + 1, Af^f 2 are the BipoSH coefficients of 
the masked sky map, correspond to the BipoSH co¬ 

efficients of the unmasked sky, W™ are the BipoSH coef¬ 
ficient of the mask itself, are the Clebsch-Gordon 

coefficients, and the term { } in Eq. (D.2) is the 9j— symbol. 
This quantifies the coupling between the BipoSH coeffi¬ 
cients of the CMB sky map and those of the mask itself. 

The underlying CMB sky may have deviations from sta¬ 
tistical isotropy, as discussed in Sect. 6.4, due either to a 
dipole modulation (L = 1) of unknown origin, or to Doppler 
boosting (L = 1) of the temperature field. The BipoSH co¬ 
efficients of such statistical isotropy-violating fields can be 
given by 

(D.3) 

Here corresponds to the BipoSH coefficients of the 

unknown but statistically isotropic CMB field. This cou¬ 
ples with BipoSH coefficients of the mask to introduce a 
mean field linear bias ask, which is estimated from 

simulations and subtracted from the BipoSH coefficients 
obtained from the masked sky. The </>lm are the spheri¬ 
cal harmonic coefficients of the field that breaks statistical 
isotropy, and Gf li2 is the shape function. Shape functions 
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for dipole modulation and Doppler boosting are given in 
Eqs. (54) and (56), respectively. 

Due to symmetries of the mask, which is largely de¬ 
fined by foreground residuals towards the Galactic plane, 
the dominant BipoSH modes of the mask correspond to 
J = {0,2},^ = 0. Hence, for all practical purposes, sig¬ 
nal is retained in the L = 1 mode itself, although masking 
modifies the shape function, now defined as the modified 
shape funtion in the rest of the text. A weighted variance 
modified shape function is defined as 



(D.4) 


where = A\^ - (-4^) mask and the weights are cho¬ 



sen such that 22, 


Here is the MSF, which can be evaluated as 




(D.5) 


The weights are then given by 


-l 
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(D.6) 
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